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ABSTRACT 


The basic outline of a digital computer (CRC-102A) program for the 
numerical solution of the approximation problem in the time domain is 
presented. The method used, in effect, is a real time convolution and 
Dirichlet series representation of the Laplace transform to effect a 
time-to-frequency transformation. The major portions of the digital 
computer program and the complete flow chart to solve the problem are 
included, 

The composition of a basic library of programs for the solution of 
all network synthesis problems and the integration of this program into 
such a library is discussed. 

The program for the solution of the roots of an E degree poly- 
nominal and the Dirichlet series expansion program were completed during 
the author's industrial tour at the National Cash Register Company, 
Electronics Division, Hawthorne, California. I take this opportunity to 
thank Mr. Henry Kent and Mr. Philip Sunday of that organization for their 
help in the writing and the debugging of the programs completed during 
my industrial tour, 

Special thanks are also due to Professor E. J. Stewart for his 
recommendations and constructive criticism with regard to the programming 
of the overall problem, and to Professors M. L. Cotton and J. B. Turner, 
Jr., for their guidance and help throughout the preparation of this 


paper. 
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l. Introduction 

The problem of network synthesis is a difficult and complex one. 
Consider the synthesis process to be divided into the two separate sub- 
problems of approximation and realization. There are many techniques 
for the attainment of the first goal, Z(p), the transfer SIRE 
Similarly, having attained Z(p), there are many procedures available to 
realize Z(p) in any one of several different circuit configurations. 
From the point of view of the electronics research engineer, a tool which 
would rapidly and accurately compute the circuit components of a network 
to provide a specified response, at a specified impedance level, in a 
specified circuit configuration would free the researcher from the tedium 
of designing the network and also from the temptation of accepting a 
"second-best" solution. The designer would feel no compulsion to term- 
inate a problem rather than face another "messy" synthesis problem. 
The modern high-speed digital computer offers the engineer such a design 
tool, if an appropriate library of synthesis programs is available. 

Possibly the greatest advantage the computer avails its users is 
that once a certain type of problem has been solved, the method of solu- 
tion may be permanently recorded on punched cards or on magnetic tape or 
paper tape. The engineer now need only recall that such a problem has 
been programmed and refer to the catalogue of computer programs rather 
than to a textbook to refresh his mind on how to solve the problem by 
paper and pencil. 

"o a comprehensive summary of approximation methods, the reader 
is referred to Trans, IRE, Prof Grp on Circuit Theory, Vol. CT-1, No. 3, 
of Sep 1954. In addition to the summary, the first paper presented by 


Stanley Winkler, has a definitive bibliography of 240 items on network 
synthesis. 





It is the purpose of this paper to present a program for the solu- 
tion of the approximation problem in the time domain and to show how this 
program may be integrated into a library of programs which use convention- 


al methods for the solution of the more common network synthesis problems, 





2. Statement of the Problem 

The approximation problem as considered in this paper is felt to be 
defined as: 

Given the excitation available or expected at the input terminals 
of a passive four-terminal network, and the response or output required, 
or the impulse response, where these parameters are expressed in graphical, 
or general analytic form; find a rational function expressed as two poly- 
nominals, where one is understood to be the numerator and the other the 
denominator, which describes the transfer function within the limitations 
on the allowable deviation of the output or on the number of elements 
to be used (but not both), specified. Provisions are to be included for 
compensation of dissipative effects. The polynominals will be expressed 
either as polynominals or as the quadratic factors of polynominals. The 


function is to be physically realizable. 





3. The Approximation Problem in the Time Domain 

Heretofore, approximation problems in the time domain have involved 
lengthy and inaccurate frequency domain expansions which do not always 
lead to realizable system functions. These methods involve complicated 


Eo At any rate, 


integrations, or at best, Taylor series expansions. 
they are not suitable for the flexible program requirements we desire. 
What we would like to have is a method which would rapidly transform the 
input-output characteristics of one network into a realizable system 
function within any desired degree of accuracy, or alternatively, an 
optimum realization for a specified number of elements. The numerical 
method of Ba ins» is tailor-made for these requirements. Furthermore, 
not only will it operate on input-output relationships, but also there 
is no reason why an impulse response which has been obtained analyt- 
ically cannot be plotted and introduced into the program at the appropri- 
ate point to produce the transfer function for all types of networks, 
i.e., predictor networks, smoothing filters, interpolation filters, com- 
pensation networks and so on. In fact, we may say if a physically realiz- 
able h(t) can be defined, we can find the system function. The designer 
must have a plot or sequence of ordinates which describe the h(t), or 
the f(t) and £ (t), and use this plot or sequence to enter hís problem 
into the approximation program. 
Briefly stated, Ba Hli's method is to find, by a process of synthetic 
division, a sequence ( Cnt A? which represents the impulse response of 
the network which in turn relates the input waveform to the desired output; 
Freddy Ba Hli, Network Synthesis by Impulse Response for Specified 


Input and Output in the Time Domain, Tech. Rpt. No. 261, MIT Res. Lab. 
of Electronics, July 31, 1953 (hereinafter abbreviated as NSIR-BH). 
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and then to calculate the psuedo-system function H*(s), from n] A by a 

simple power series expansion. The H*(s) is a polynominal in s which is 
the ratio of the En whieh is normally considered to be the expression 
of the system function. In effect, we convert the convolution integral 

into its component operations and use an iterative substitution method. 

We shall not reproduce the theoretical reasoning and proof which was so 

neatly done by Ba Hli, but we shall demonstrate the ease and simplicity 

of the method with a numerical example, 

Before proceeding with a complete example, let us show how the 


synthetic division process yields the impulse response in a couple of 


trivial cases, the results of which are well-known. Consider Figure l, 


Unit Step 


ít, c] りー T: 


Unit Ramp 


fe ce} sU O 5, L0, a 





A Trivial Case 


Figure 1 








the input f(t) we take as the unit step whose Laplace transform is 1; 

f (t), the output, is the unit ramp. Its transform, of course, is 2. 
If we sample the input and output at intervals of .2 of a time unit, the 
sampled ordinates for f (t) are uU UE and for £ (t) 


[+2,.4,.6,.8,1.0,1.4,1.6 ET ie Now dividing f (t)) in the manner 
EY 
DE Ls AS eee 
sllustrarad: 1,1,1,1,1,1..../.2,,.6,.8,1.0,1.2, 1.4.02... 





we have UN m ,.2,.2,.2,.2,.2,..... and if we divide Int A by At, we pro- 
duce the sequence 1,1,1,1,1,....again the unit step. This, of course, is 


completely analogous to the process in Laplace transform operations, viz: 


Lee] Lhe) = Lhe 


2 


We can apply the process in reverse, for example, (£, (09) X {n} A 
= (£,(6)), if we use the unit doublet as inj a we should differentiate 
the output. Consider the sequence below which defines a sine wave, with 
Acre. ol 
{f} = .1 .199  .296  .389  .479  .564  .644  .717 


{h}, = 10  -10 


1 1.99 2.96 3.89 4.79 5.64 6,44 7,17 





21 -1.99 -2.96 -3.89 -4,79 5.64 6.44 
if 1 .99 .97 .93 .90 .83 .80 .73 


i£ (Cont'd) .783  .841  .891  .932  .964  .985  .997 


7.83 8.41 8.91 9.32 9.64 9.85 9.97 


7.17 7.83 8.41 8.91 9.32 9.64 9.85 


\f (Cont'd) 66 .58 | .50 . .41 .32 22008 12 








and O) is a sequence which defines the cosine which is as expected 
since the unit doublet acts as a differentiator. Thus we demonstrate 
the ability of this (ni A to perform as an operator equivalent to H(s) 
in the frequency domain on functions defined as sequences in the tíme 
domain. The transformation of fh} A to H*(s) follows logically from 


the definition of the convolution integral: 


t 
f (t) = E, f (t) h (t -T) aT 


Ba Hli 2 has shown that H*(s) can be represented by the Dirichlet 


power series, 
o 


H* (8) = 7 a exp (-T 8) 


nzl 
where the a, are the areas in our Í h} A sequence, and the T. are de- 


fined by 


E az nåt - AE 
n 2 


> 
that ís, the | are measured from the origin to the center of the sampl- 
ing interval under consideration at the moment, 

We may approximate H*(s) as closely as we please by taking our 
sampling points at smaller and smaller intervals. 

To demonstrate the complete process and provide a basis for diseuss- 
ing some of the more obscure points, let us take another example and 
follow it through more closely. Figure 2 shows the £,(t), the desired 
£ (t) and the resulting h(t) required, (At # 0.1): 

(Et = 1, PNEU ti. I. Gi fo. O, 
(£ (E) 1, SEE A EO 
WE, 1.1041, 21, 0, 0, 0, 0, 


In tabular form, we have 





SEE 








zu 





EE 


coh 





rt 
| 

= ++ 
! 


EE 





E 5 E seit li. 
- de ET — + 
E ゴー E + — H 
E u 
[all 


He Dc 
ELA 
は H-- hid 


HRP PP 








.05 1 
LO 1 
25 el 
235 el 
.45 XT 


Substituting into our expression for H*(s), we have 


2 3 4 O 
Hs)»  .l[r-.o5, GHH- mer , C93). . |] 
(.15a)^ (.15s)" (sn: 
+ œl E - „l5s + 21 - 3! * à? - T" 
(„258% ( 258)" (,25s )* 
+ aL 1 = .255 + 21! = = 3! + 4! T - 
( TP ( e (,35s )* 
* en f x 。35s ~ = 21! = 3! + 4! d Be 
2 3 4 
film ne GE. A, use LL] 


= .5 - ‚1258 + .020625s* - .00255208s? 4 .0002517968s* - ....... 


If we compare this result with that given by Laplace transform methods, 


we have 
3 t, ce) 
Jr (el à =, L1 - exp (-1/28)] * 


A [nt )] Ay, (で) E - exp (-1/28) 


u 
の コー 


| 1 - exp (-1/2s) | 


and the power series expansion for \(h(t)| is 
.5 = .1258 + .020833s^ - .0026041665” + ‚00026051665 
which shows the first two terms by Ba Hli's method agree exactly, a 1% 


error in the third term, 2% in the fourth, and about 4% in the fifth. 
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Before proceeding let us see how this accuracy may be improved. 
Clearly increasing the number of terms taken and maintaining the same 
interval will have no effect on the formation of the second, third, and 
fourth degree terms whose accuracy we wish to improve and intutively one 
feels that increasing the number of sample points will; therefore let us 
increase the number of sample points to, say, 8. 


With f (t), £ (t) and h(t) as before we tabulate, for At = .0625: 


Tn ?n 

03125 .0625 
.09375 .0625 
.15625 .0625 
.21875 .0625 
.28125 0625 
34375 .0625 
40625 .0625 
46875 0625 


2 3 4 
H(s) = .0625 [1 - .031258 Bel . (031238). , (0931238). . ,, | 


2 3 4 
062s [1 = .093758 + £09375)" _ 14093758)" _ 1.093736)". > 


and so on; the resulting expression is: 

H¥(s) = .5 - .1258 4 .020635719s^ - .002631105s? «4 .000257032861s^ - ... 
and the new expansion for 8 sample points has produced errors of 0% in 
the first two terms, .95% in the third term, 1% in the fourth and 1.3% 
in the fifth. Therefore, the increase in the number of points has 
obviously resulted in an overall improvement. A further increase to 

10 points, produces accuracies of 0%, 0%, .25%, .5%, .8% and 1.25%. 


The foregoing results are tabulated in Table 1 for the purpose of 
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easy comparison, 

We recognize that the errors in the higher order terms are caused 
by replacing the Laplace-Stieljes inter pal” by a truncated series summa- 
tion. This is equivalent to saying that we have not computed the true 
"moments" of a, about the origin. 

Consider Figure 3, we have here a triangular impulse response. Let 
us calculate the second and third moments of this curve and compare the 
results with the power series expansion of the Laplace transform and 
also the Dirichlet series expansion. We let b, be the second "moment", 


(it would be more precise to call this by a different term since it is 


actually E. , Where M, is the moment), b, equals the third; thus, 





1! : i 
1.0 
Es h(t) =t, 0z t «€ 1/2 
l-t, 1/2£t<1 
0, elsewhere 
h(t) $2 
0 ~ 
0 の 1 E 
A Triangular Impulse Response 
Figure 3 
1 (3) 


NS IR-BH ? Pg. 15. 
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LSZ°1- Zct8'- 40065” - ROSZ- 40 20 | 102243 3U9) 194 
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1/2 1 
b, = ーー | E E + » (1-x) x^ax| 


1/2 1 





1 
4 
0 3- 1/2 zi al 


=  .03645833 


1/2 1 
b, s L| 2 EM n (1-x) x dx | 











00781250 


The trancendental Laplace transform for this function is: 


H(s) = iy [1 - exp ers]: 


whose power series expansion is 


H(s) = .250 - 0.1250s «4 .03645833s^ - .00781250s" 


The Dirichlet series results from 


Th a, 
x . 02 
a .06 
ae . 09 
3 . 06 
9 . 02 

-T s 


Ee a, € dl 
n=l] 
Therefore, considering only the b, and b3 terms, we compute: 
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2 2 2 2 2 
er 02 x —— 


2 2 2: 2: 2. 25 
a > (.0002 + .0054 + .0225 + .0294 + .0162) 
= .03685 
Similarly, b, is computed as 
bz = .008008 
Therefore: 
"Moment" Laplace Ba Hli Approx. 
b, 03645833 03645830 .036850... 
b, .00781250 .00781250 .008008... 


The obvious method of reducing the error is to take a smaller in- 
terval and thereby approach the true integration to a closer degree and 
since this will be a one time only calculation, we are not adamant about 
making extra computations to improve the accuracy. We also notice that 
by making our intervals smaller we reduce those errors which may arise 
from the process in which we obtain the {ht AS For a more complete dis- 
cussion of the reduction of the inherent error in Ba Hli's method see 
Appendix II. 

Having obtained a psuedo-system function, it yet remains to deter- 
mine its equivalent as the quotient of two polynominals, the degree of 
each as yet unknown. We equate our psuedo-function, H*(s) to an express- 


ion of the form 


m 
P(s) Po + P$ + ... Ps 
n-1 n 


do + Are + cee Neg? * 19 


H*(s) = 


and solve for the unknown p's and q's, since these coefficients express 
a ratio we may arbitrarily set any one of them we choose, (normally q,) 
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equal to one without any loss of generality. 
We notice that the series of equations which results from this opera- 
tion will always have the form: 


AA oer 


Po = hodg + ha, * hgg, 


0 s hado + h,q, + hid, + hod, 

0 = b dg + hq] + DP - hiq, + hod, 
ner Fd 6 C838 
ME UE. up. * 545 + BT 159, 


For m = 2, n =4; r, the relative degree between P(s) and Q(s), is 
obviously two. We see that we have seven equations in eight unknowns, 
which is as it should be for the expression of a ratio, and hence, re- 
quiring that we set dg (or any other coefficient) equal to some conveni- 
ent value. The solution of this portion of the problem is considered in 
detail in Appendix III. 

We choose r, the relative degree of numerator and denominator by 
observing the behavior of h(t) around 0. Since we know from the initial 
value theorem of Laplace Transform Neary: 

lim h(t) = sĦ(s) 

t0 s ——> »o 
we may say that the function must behave as 


Pa? 


S 





as t approaches zero, 


In. F. Gardner and J. L. Haines o). page 265. 
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Thus, if h(t) starts out as a step 


P(s) > l/s as s approaches infinity; or is h(t) starts 


out as a ramp P A Wee as s approaches infinity. 


To program the computer to sense this limit and thus decide whether h(t) 
is an impulse, step, ramp, or higher order function, we first determine 


ha, i.e., h(0). Since in some cases this may result in division by zero, 


0’ 


we perform a first-order backward interpolation, using hy and h, to find 


h This is approximate but certainly is sufficiently accurate to deter- 


0° 
mine whether or not h(t) starts as a jump function or as a ramp, which is 
the primary decision in the "selection of relative degree routine” (see 
Appendix I, page I-1). If it is decided at this point that h(t) is a 
first-order approximation to a ramp function, we must further decide 
whether it is truly a ramp or a higher order function, i.e., War, 17: 
etc. These decisions are based on the fact that the nt derivative of an 
n-order polynominal is a constant, the process is flow charted for 

*l — r £ -4 on page I-2-1 of Appendix I. The routine may be easily ex- 
tended to encompass all possibilities which may arise. 

We know that we can obtain any desired degree of accuracy in the 
terms of our power series expansion and further that we can correctly 
determine the relative degree between our numerator and denominator poly- 
nominals., We have shown that we can compute the coefficients’ of our H*(s) 
and we know that since this is a ratio with arbitrary coefficients we may 
multiply the expression by the impedance level factor to set any required 
impedance level. But what of the approximation error, i.e., the error 


caused by not taking an infinite number of terms to approximate tran- 


cendental functions? 


tpa Hli has shown that in the case where the Laplace transform of 


h(t) is rational, it will be recovered exactly. NSIR-BH,(3) pgs. 41-44. 
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Ne 5 -1 + 1.5(5) + 2 = 13.5, which we take as 14. 
Therefore, we can write 


S s? + + 51 + s13 
un eFt P, PI P13 11 27 


15 = ho + his e soo h, 18 + see h, 78 


dy + 9,8 + ans" +... + TO * q448 
this expression leads to 23 equations in 24 unknowns, one of which, custo- 
marily dg» we set equal to 1, without error. To avoid filling sheets of 
paper with useless symbology, we shall terminate this example here and 
turn to another pair of functions. 

We shall compute, using the program, approximate impedance functions 
for a pair of input-output time funetions and check the computer solution 
analytically. To do this, we choose the functions illustrated in Figure 
4, 

£,(t) is a rectangular pulse 

f (t) = 1 0o£t 三 1 
0 elsewhere 


£ (t) is an inverted sectioned parabola 


f (t) = oe 0=r*ı 
= ee ler *2 
E 2£t #3 
= 0 elsewhere 


for At = .l, the following sequences obtain 

O] = 1, l, l, l, l, l, l, l, l, l, 

{f (t) = .005, .020, .045, .080, .125, .180, .245, .320, .405, .500, 
.390, .660, .710, .740, .750, .740, .710, .660, .590, .500, 
.405, .320, „245, .180, .125, .080, .045, .020, .005, .000, 


Synthetic division of EO yields: 
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In} a = ・005。 .015, .025, .035, .045, .055, .065, .075, .085, .095, 
.095, .085, .075, .065, .055, .045, .035, .025, .015, .005 
which describes h(t) as the triangular function also shown in Figure 4. 


We know that the Laplace transform from ft) is: 


Alf] = 3 a-e®) 


and the Laplace transform for f (t) ms 


S 


0 [=] 


- =s 2 
Thus H(8), which should be EXO] / IE) x pe en 
the transform of a triangular ae 


We expect H*(s) as computed by our program to be a good approximation 


-s— 2 

to the power series expansion of LÀ < 
NH e 7 2 2 
ュー を | = 1.000000000 - 1.000000000s + .583333333s 
- .250000000s? - .086111111s* - .025000000s? 


6 


+ .006299603s .001505423s^ +.... 


The program approximation with single interpolation, i.e., At/2, produced: 


Rx(s) = .999999999 - 1.000000975s  - .583340832s^ - .250007019s” + 


.0861150773* = .0250016708” 5 .0063001725° = .001405587s/ +. 


If we now solve for the impedance function using n s 3, and r = -2 


(since h(t) is a ramp function in the vicinity of zero); 


0.999999998 - 0.049971105 
z(p) = 0.99999999 s 


or; 
1.0 + 0.950029868s + 0.3666899615s^ 4 0.062506119s? 


x» (-,049971105) (s - 20,011564640) 
(s + 2.277182620) f(s + 1.79464120) + 1.950589991*] 


er and kaes 00) Appendix A, No. 6.08. 
2 (6) 
Gardner and Barnes , Appendix A, No. 6.07. 
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The inverse Laplace transform of this function is: 


ana 7 e °4/718202t „465570186 e «'/2*495120€ .in (1.95058999t-81.217?) 


The plot of this inverse superimposed on the desired h(t) is shown in Figure 
5a. 
If we desire to improve this approximation we go to a 4-pole network. 


The program computes the system function as: 


2674145558" - 8,16227058 + 135.9105227 


Dn + 10.598284918” + 51.1404969s* + 127.7483847s + 135.9105230 


which we expand as: 


10681279249 1 s + 23,98432161 - s + 20.00984647 N 


| [(s43.017294581)*«(.932721456)^] [(s+2.281847873)%+(2.901655110) | 
and determine the inverse Laplace transform to be: 


-3.017294581t 


24.010971222 e sin (.932721456t + 2.547”) 


a o SIC A (2.901655110e4 9.295°) 


This inverse is plotted in Figure 5b. Note the improvement with the 
‚addition of the other pole. The approximation can, of course, be further 
improved by adding still more poles. 

Suffice it here to say, that the existing program will solve for the 
impedance function, given the time-series data describing the input and 
output waveforms. 

Before concluding this section, we estimate the required solution 
time for the total approximation problem in the time domain, 

We assume the following conditions obtain as a basis for our estimate: 

l. h(t) is not a closed function, i.e., K = 4 (see Appendix II). If 
h(t) is closed, that is, if h(t) is completely EL in 4 or less time 


units the running time for the synthetic division routine and for the com- 
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putation of the psuedo-system function, H*(s), will be proportionately 
reduced, 

2. H(s) starts as a step function, i.e., r = -l, This is the worst 
case we encounter as regards computation time for the psuedo-system funct- 
ion, since we must compute 2(n + l) + r terms of the expansion of H*(s) 
in order to solve the system of simultaneous linear equations. 

3. The interpolation decision is such that one intermediate inter- 
polation calculation is performed. This is the intermediate condition, 
and requires twice the computation time for a no-interpolation decision 
and one-half the time for the four-point interpolation. 

We consider the overall program (for time domain approximation) to 
be sub-divided into three, (four if Z(p) is desired in factored form or 


if compensation is to be provided) major computations. Thus we have: 


Sub-division CRC time Word-times 

* 
Synthetic Division to 5 min* 72 x 10^ 
produce {hl A 
Dirichlet power series n x 3 min (n-1) x 44 x 10* 
expansion H*(s) 

3 3 4 
Matrix Inversion to -064 n min 96 n x 10 
compute Z'p 
Factor Z'p to apply (2n-1) x 9.4 min (2n-1) x 144 x TN 


loss compensation 
*Estimated, all other values measured, 
For a 10-pole network, this gives a solution time of 278 minutes; for 


a l5-pole, 538 minutes; for a 20-pole, 945 minutes. 
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4. The Network Synthesis Library 

Figure 6 shows the block diagram of a network synthesis library for 
the automatic computation of the approximation functions for a variety of 
filters. This diagram is by no means considered to be all-inclusive, it 
is only intended to show how such a library would be assembled. 

It should be obvious from the diagram how our various blocks are used 
together to provide a large number of filter characteristics from a rela- 
tively small number of programs. All blocks would have common input, 
output and carry* locations, and programs represented by the blocks would 
have bootstrapping links to call in the next block in order to proceed 
with the computation. 

In addition to the "frequency" filters shown, we would consider it 
desirable to provide in the basic library: 

a. Tschebycheff Stop-Band Ripple Filters 

b. Tschebycheff Equi-Ripple Filters 

With regard to the latter, Byron J. AS has developed an inter- 
esting iterative technique which seems to be much simpler than Alexander 
J. Een method for avoiding the elliptic function. 

Note the provision for the realizibility check. It is considered that 
the inclusion of such a routine is mandatory in order to save computer 
time. This routine would compute E the number of elements required to 


*We use the term carry to characterize that data which is not used 
in the present computation but must be "carried" for use in later routines. 


“Byron J. Bennett, A Note on Circuit Synthesis, p. 61,616) 


a. J. Grosman, Synthesis of Tschebycheff Parameter Symmetrical 
Filters, Proc. IRE, Vol. 45, No. 4, pP. 454-474.(13) 


(13) (4) 


ne Grosman's equation (33) or see Darlington : 
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provide the required selectivity parameter (or discrimination parameters) 
and meet the ripple requirements, if n is specified and the requirements 
demand more elements than i allowed the computer should so indicate to 
the operator. It is estimated that this information would be available 
within not more than ten minutes from the start of computations, and 
therefore establish a doctrine that if such a program is used, the opera- 
tor must wait until this point in the program is passed before allowing 
the computer to proceed unattended, and by so doing be assured that the 
results will be realizable as a physical network, and that computation 


time will not be wasted. 
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5. Conclusions 

The results of the tests on the routines programmed to date indicate 
that a library of network synthesis routines employing conventional tech- 
niques is practicable. 

The use of such a library should result in considerable savings in 
time and effort for those people concerned with network design. 

The accuracy obtainable is more than sufficient to guarantee excellent 
approximations (and undoubtedly synthesis) if the number of elements to 
be used is not severely restricted. 
| The time required for solution is not excessive in comparison with 
the magnitude of the problem and the accuracy of the results. The employ- 
ment of the library of programs will require a working knowledge of the 
computer, unless magnetic tapes are made up from the “building blocks" 
to provide specific types of network designs. In this case, detailed 
instructions may suffice to permit persons not familiar with the computer 
to use the programs to good advantage. 

The installation of the CDC - MK 1 Computer will enhance the value 
of this library despite the convention problem. The increased speed of 
computation and shorter programs which will be possible due to the "built- 


in" floating point arithmetic will make such a library even more valuable. 
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APPENDIX I 


It is the intention of this appendix to explain with a minimum of 
theoretical reasoning the overall flow chart of Ba Hli's method, 

The following inputs are required: 

I. ít (t) - time sequence of output function (in octal or decimal). 


25 i£, (t) - time sequence of input function (in octal or decimal). 


3; N - number of discontinuities. 

4, N. - number of points of inflection. 
5. N - number of elements. 

6. t - number of f(t) ordinates used, 


7. Octal or Decimal input data indicator. 

Annex 1 shows a block diagram flow chart which considers only inputs 
and outputs to major program subdivisions, and major logical decisions 
required, 

We allow for either decimal or octal digit input sequences, in either 
case they are converted to binary floating point. If the number of ele- 
ments (N) is not specified, compute AS = 1.5 Ny + N. (see page 17,main 
text). 

The synthetic division process is self-explanatory, the quantity 
"t", the number of f (t) ordinates, is required in this block. 

Annex 2 shows a detailed flow chart for the determination of "r", 
the relative degree between numerator and denominator. The first order 


approximation used to compute h, is considered to be sufficiently accurate 


0 
to ascertain whether or not h(t) starts from a finite value, or from 
zero. The "O" used in the decision boxes is actually slightly greater 


than zero, about .025, to avoid any erroneous decision that ho is greater 





than zero, when it is only greater due to the error introduced by the 
first order approximation. Note that we have deviated from Ba Hli's 
method of determining the behavior of h(t) near t = 0 in that we compute 
derivatives. 


We can now determine, D the degree of the denominator, and N 


d’ d’ 
the degree of the numerator of Z(p) and the number of terms of H*(s), 
Hye we must compute to solve for Z(p). 

We can make a somewhat arbitrary decision as to how accurately we 
should compute our H(s), i.e., what subdivisions we should use in our 
"moment" computing scheme (see Appendix II), depending on the number of 
terms we shall have in our denominator. 

The next two appendices give the program details for the solution 
of H*(s), and Z'(p). 

The solution of Z(p) is now complete, unless we desire to compensate 
for dissipative losses. This is easily accomplished after factoring the 
numerator and denominator polynominals. We preshift our pole and zero 
pattern to the right an amount "x", unless this shift causes a pole or 
zero to move into the right half plane and thereby violate the conditions 
of physical realizability. We compute 

X = R/2L 

2L/R = 1/« 

wL/R = の /2 へ ん 

2Q /o - A 

N a Q/rf 
and thus we compute «a, as a function of Q and f, where f is the lowest sig- 


nificant frequency involved and is entered by the designer as input data. 
















Determine r 
5 +r+ AS ee D 


D, -r -N 


d 


D + Na 


Interpolation Decision 


Hx(s) 


See Appendix II 





Convert to BFP 
Determine AS if N =0 


Synthetic Division 
to produce Üp 


See page I-2-1 





Pm(p) の 
Qu) m 


See Appendix III 


If Q entered, 


determine A 





Factor P_(P) 






Preshift by A 








Factor Q (a) 
Preshift by A 


Realization Programs 


OVERALL FLOW CHART 
OF 
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APPENDIX II 


CALCULATION OF H*(s) FROM THE Ih! SEQUENCE 


A 


We are now required to form the power series expansion of: 


As was mentioned in the main body of this paper, replacing the Laplace- 
Stieljes integral by the summation indicated above, resulted in an increas- 
ing error in those terms higher than the Iter degree. We recognize that 
this error results from the failure to compute the exact moment of bn} a 
about the origin. 


The true qa moment of an area defined by f(x) is: 


50 
Ma = E x f(x) dx 


which we have replaced by the summation 
S ZUR 
M Capp) = 2 E A. 
iel 
where A, = incremental area or na, and S's NK (K= 1,2,3.... ) 

For our program, we limit K to 4; that is, we require the impulse 
response to be defined over the range 0 to 4 time units. The input and 
output functions may be sealed so this will be adequate for most problems. 

It might be well before proceeding to take a close look at the manner 
in which we form the coefficients for the Dirichlet power series expansion. 

We may proceed in such a manner as to form one term at a time, that 
is, referring to Figure II-l, we may sum first down the j = 0 column, then 
the j = 1 column, etc. The operation count for this procedure is: 


N "a" for the first term, 


LL tl 





N "a" 4 N'"m" + "m" for the second (if we factor out the 1/2N until 
the sum is complete), 


N "a" +4 N(j + 1)'m" + (j + 1)'m” for the nn term. 


a total of 
j+l 
(j + 1) Wa" + Nel 2 (i+ 1) "m" 
isl 
(j + 1) fn LE * (= + 1) "m" | +( He + 1) m") 


where "a" stands for floating point addition and "m" stands for floating 
point multiplication, 

For N a 40, j = 20; this becomes 840 "a" + 11040 "," and if we assume 
that floating point addition requires twice the time that floating point 
multiplication does, a equivalent total of 12,720 floating point multi- 
plications. 

If, however, one proceeds out the rows and forms all the partial sums 


associated with A then A,» etc., for i = l1, j = 0, to j, the count per 


1? 
row is: 

Tae 2 me Va "a" + "m", at + "m", .... "a" + "m" and the 2N in 
the denominator is included. 

The count is the same for all rows, therefore, the total is; 

N(3+1)("a''+'"m'"), which for N = 40, j = 20, is 840 "a" + 840 "n", or 
based on our previous assumption, a total of 

2,520 floating point operations. 
There is, however, some additional "call-up" and "put-away" address compu- 
tation which degrades this efficiency somewhat. Nonetheless, it seems 
quite evident that a saving of at least 50% can be effected in this way. 


With this saving we can afford to reduce the error in the final co- 


efficients by introducing an interpolation routine which effectively 


11 - 2 





doubles or even quadruples the number of sample points used. 

The interpolation precedure makes a straight line approximation be- 
tween sample points, and computes the value the intermediate sample points 
would have. Therefore, it is able to use a smaller increment on the 
abscissa and thus approaches a true integration more closely than would 
be otherwise possible. The results of some runs using various "degrees" 
of interpolation are presented in Table II-1 for comparison purposes. 


The flow chart and program are appended hereto. 
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The factorials which appear in H(s) are yet to be supplied, 
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ADD 


0300 
0301 
0302 
0303 
0304 
0305 
0306 
0307 
0310 
0311 
0312 
0313 
0314 
0315 
0316 
0317 
0320 
0321 
0322 
0323 
0324 
0325 
0326 
0327 
0330 
033] 
0332 
0555 
0334 
0555 
03 36 
03537 
0340 
0341 
0342 
0343 
0344 
0345 
0346 
0347 
0350 
0351 
0352 
0353 
0354 
0355 
0356 
0357 


INST 


M1 


3000 
3000 
2010 
0556 
3000 
2020 
0541 
3000 
0542 
0536 
0510 
0555 
0120 
0555 
0523 
0555 
3000 
0545 
0572 
0523 
0564 
3000 
0523 
0524 
3000 
2000 
2000 
0523 
0521 
3000 
2000 
2000 
0523 
0521 
3000 
2000 
2000 
0523 
0534 
3000 
0572 
0525 
0561 
0525 
0564 
3000 
0525 
0524 


2100 
0532 
2001 
0564 
0562 
2100 
0532 
2001 
0564 
2100 
2100 
0533 
0534 
2100 
0532 
2100 
2100 
0564 
0565 


M? 


0550 
2000 
0505 
0545 
0311 
0310 
05h 5 
0311 
0545 
0555 
0314 
0551 
0523 
0321 
0530 
0547 
0413 
0350 
0326 
0521 
0562 
0333 
2000 
2002 
1516 
2000 
0521 
2000 
2002 
1516 
2000 
0530 
2000 
2002 
1540 
2000 
0527 
0524 
0547 
0413 
0356 
0521 
0562 
0522 
0563 
0365 
2000 
2002 


DIRICHLET POWER SERIES EXPANSION FOR PSEUDO-SYSTEM FUNCTION 


REMARKS 


Load buffer 
and clear putaways 
If TS 1 up, 0305 (1) 
No, Set t to 1 
and skip to 0311 
If TS 2 up, 0310 
No, Set t to 2 
and skip to 0311 
Set t to 5 
N exp + t to N’ exp 
Set hy pickup 
に to SA 
If t > 1, 0321 
TY to TA 
1 toH 
Skip to 0413 
If t > 2, 0350 
No, if SA > "0", 0326 
No, TF/2 to A 
Skíp to 0333 
Load TF 
Load TE 
TF - TE 
(TF - TE)/% 
(TF - TE)/4 to A 
Load TF 
Load à 
TF - A 
(TF -A)/2 
to TA 
Load TF 
Load A 
TF 44 
(TF + A)/2 
to TB 
TF to TE 
2 to H 
Skip to 0413 
Is "Sa" > "9" 
No, TF/h 
Los 
TF/2 
tod 
Skip to 0365 
Load TF 
Load TE 


(1) The test switch function would be performed by the interpolation 
selection in the overall synthesis program. 
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ADD 


0360 
0361 
0362 
0363 
03 64 
0365 
0366 
0367 
0370 
0371 
0372 
0373 
0374 
0375 
0376 
0377 
0400 
0401 
o402 
0403 
040% 
0405 
0406 
0407 
0510 
0h11 
0412 
0413 
0414 
0415 
0416 
0417 
0420 
0h21 
0h22 
0423 
0424 
0425 
0426 
0427 
0430 
0431 
0432 
0433 
0434 
0435 
0436 
0437 


M2 


3000 
2000 
2000 
2000 
2000 
0523 
0522 
3000 
0521 
3000 
2000 
2001 


3000 
2000 
2001 
0523 
0521 
3000 
2000 
2001 
0522 
3000 
2000 
2001 
0523 
0543 
0531 
0534 
3000 
2000 
0535 
3000 
2000 


0520 
2100 


051 


0550 
0520 
3000 
2000 
0600 
0530 
3000 
0433 
2006 


M? 


2100 
0534 
2001 
0532 
2001 
0564 
0563 
2100 
0562 
2100 
05 34 
2100 
0563 
2100 
0534 
2100 
0564 
0562 
2100 
0534 
2100 
0563 
2100 
0534 
2100 
0564 
2100 
0572 
0575 
2100 
2001 
0576 
2100 
2001 
0516 
2100 
2100 
2100 
2100 
0571 
0561 
2100 
2001 
0601 
0571 
2100 
1505 
1677 


MS 


1516 
2000 
0522 
2000 
0521 
2000 
2002 
1516 
2002 
1516 
0530 
0571 
2002 
1540 
0527 
0570 
2000 
2002 
1540 
0526 
0567 
2002 
1540 
0525 
0566 
0524 
0547 
2000 
2002 
1540 
0531 
2004 
1507 
0520 
0433 
0457 
0544 
0427 
0433 
2000 
2004 
1500 
0530 
2000 
2002 
1540 
2006 
0442 


11 = 1-2 


DIRICHLET POWER SERIES EXPANSION FOR PSEUDO-SYSTEM FUNCTION 


REMARKS 


TF TE 
(TF - TE)/% 
to A 
(TF - TE)/8 
to A 
Load TF 
Load A 
TF <A 
Load A 
TF 24 =- ô 
(TF -4 -A)/h 
to TA 
Load A 
TF- À 
(TF - A)/h 
to TB 
Load TF 
Load A 
TF +A 
(TF 4 A)/h 
to TC 
Load ^ 
TF + ム + み 
(TF + ム A+ み )/4 
to TD 
TF to TE 
Set H to 4 
Load SA 
Load Mo 
SA + 2 
to SA 
Load N° 
Divide 
Result to M 
Reset putaway 
If M30, 0457 
No, clear j tally 
If j 0, skip to 0427 
No, skip to 0433 
Load TA 
Load M 
TA X M 
to TA 
Load summation 
Load TA 
Accumulate 
Build putaway 
Build putaway 





ADD 


OLLO 
0441 
0442 
0443 


INST 


00 


- 0536 
02 


00 


M2 M? M? 
0442 0532 2006 
2006 1677 0443 
2000 2100 0616 
2001 2100 0617 
0544 0532 0541 
0h33 0512 0433 
0540 05h 0h25 
0547 0532 0547 
0547 2100 0453 
0314 0512 0314 
3000 2100 0314 
0527 0570 0530 
0526 0567 0527 
0525 0566 0526 
3000 2100 0413 
0514 0516 046% 
0600 0601 1310 
3000 2100 1011 
0532 2100 0544 
0532 0573 0517 
0620 0621 2000 
0544 0532 2106 
2006 2100 0470 
3000 2100 0471 
0441 0537 2001 
0517 0560 2004 
3000 2100 1507 
2000 2001 1310 
3000 2100 1011 
05144 0532 O544 
0551 0544 2000 
0517 0560 2004 
3000 2100 1500 
2000 2001 0517 
0464 0512 0464 
0540 O544 046% 
0000 0000 0000 
Available 
0000 0001 0000 
2000 2005 2000 
0002 0002 0000 
3000 2100 0302 
0602 0603 0000 
0000 0000 0001 
1777 7777 0000 
Used for tempoary storage 
0000 0000 0000 

0000 0000 


0000 


Teles ss 


DIRICHLET POWER SERIES EXPANSION FOR PSEUDO-SYSTEM FUNCTION 


REMARKS 


Build putaway 
Build putaway 
Putaway 
Putaway 
Add 1 to j tally 
Modify summation pickup 
If J > 3, repeat to 0425 
No, reduce H by 1 
If H > 0, skip to 0453 
No, modify hy pickup 
Repeat loop 
Shift TB to TA 
TC to TB 
TD to TC 
and repeat loop 
Reset h, pickup 
First term pickup 
Convert and print H, (2) 
Set j tol 
Load F! 
Pickup H, term 
Extract inet bit of j 
If j odd, skip to 0470 
No, skip to 0471 
Change sign to minus 
Load F! 
Divide by F! 
Load H 
for conversion and print out! 
Increment j by 1 
Load "j" 
Load F! 
Multiply for (F 4 1)! 
Putaway 
Modify H, pickup 
Have J terms been divided by 
appropiate factorial; no, O64 
Yes, halt 


hy pickup reset 
Putaway clear routine 
Modifier for 0502 
Putaway clear routine 
Puteway reset 

Putaway clear routine 
M! and M? extractor 


Sign extractor 
J 





DIRICHLET POWER SERIES EXPANSION FOR PSEUDO-SYSTEM FUNCTION 


ADD INST M? M? M9 REMARKS 

0541 00 0000 0000 0002 2 - for H 

0542 00 0000 0000 0003 3. for H 

0543 00 0000 0000 0004 h - for H 

0544 00 0000 0000 0000 j tally 

0545 00 0000 0000 0000 t storage 

O546 00 0000 0000 0001 1 

0547 00 0000 0000 0000 H storage 

0550 26 2100 2100 0600 Putaway clear routine 
0551 00 0000 0000 0044 Constant for 0476 
0552 3, 2004 2000 2000 Putaway clear routine 
0553 00 0000 0000 0000 Available 

0554 00 2100 2100 1000 Testword for clear 


0555 - 0577 Used for temporary storage 


0600 - 1000 Used for storage of results (may be reduced to provide 
a minimm of 2(J + 1) cells. 


(2) The conversion and print routine is linked in this program for 
test purposes Only and should be replaced with a boot-strapping 
routine to call in next portion of the overall program. 
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| Set K—4 


Up 
Clear Tally 
3 > *(320) 
Initialize B test Command 
\ 
+ o : 
Yes No Nexp Vz N exp 
"nal ss SA 
prep exp 
rn — ~ 
No Yes 
Add 1 to tally 
B- 2->B 
No Yes 
N 


1>H 


9 50 » ta11yN 1 





Kel —>K 


TF- A —>TA 
2 


TF å —> TB 
2 


— > TE 
2 —>H 


K - l>K' 
Clear Putaways 
Set hy pickup 


0413 





SA + "2" ——> SA 
SA + N'-—» M 
Reset PA 





Set A; pickup 


Pickup H_ 
Convert & Print 
Set j al 
I > pP?! 







Pickup H.->2000 
Extract last bit 
of 





Yes No 


j odd? 


ex(-) into 2001 


Div by F! 





Conv & Prt 
H-1 ->H add l to j 
j —> "y" 
njx F3 -> F; 
mod A, pickup 
No Yes 
mod h, pickup 
> No Yes 
J >j 
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APPENDIX III 


MATRIX INVERSION 


As was noted on pages 14 and 15 of the main text of this paper, the 
final solution of Ba Hli's method led to a system of simultaneous linear 
equations. While it is certainly possible to program the solution of a 
system of equations by the straight forward elimination technique, this 
method is awkward when many equations are involved (and for precise 
synthesis of transcendental functions, there must be many). The computer 
must do all the address computation based on N, the number of poles (or 
the degree of the denominator) and r, the relative degree between numera- 
tor and denominator. This method for n equations in n unknowns, neglect- 
ing the address computations, would require ne divisions to eliminate the 
first unknown, a. to eliminate the second, etc. In all, to reduce 
the system to the desired "one equation ín one unknown" (we are dealing 


with square matrices), requires: 


n 
( )( ) 
/ = P met fene operations* 


izl 
If, and only if, we have saved the intermediate equations, we may 


complete the solution in 


2 
n +n 
2 





operations 


or for the complete solution, 


NO 


n 
ds: n + operations 


ola 


Such a procedure would require an inordinate amount of storage 


*The term "operation" is customarily taken to mean multiplications 
or divisions. Addition, subtractions and address modifications are not 
included in operation counts. 


Lp 





space and worse, intricate address computations. 

For simplicity of programming and best relative speed, we would 
suggest that the denominator matrix, i.e., that group of equations which 
are equated to zero, and from which we determine the q's, be solved by 
the elimination form of the mer 24 and solve the remaining equa- 
tions for Pg» Pit Pp which incidentially will always be explicit in 
terms of the q's, by simple algebra. 

We will describe this highly efficient method of computing an in- 
verse matrix and demonstrate its simplicity and economy of operation 
by means of a simple example. Briefly, by way of introduction, we have 
the situation: 

[n] [Q] = [h], if we multiply by the inverse m 

[y^ [i] R] = B] [h] we have 

Q] = a]! [n] , the desired result. 

Our immediate problem is to find a rapid means of computing the 
inverse of H. The general method of solution can be summarized: 

(1) [H I] This operation can be carried out by our 

(2) fu”! H s I| suggested technique operating solely on 

an] the unit matrix. 

We perform this transformation on I by bringing in H one colum at 
at time, multiplying by those elements which have been transformed by the 
introduction of the columns preceding. 

ea: B. Dantzig, Computational Algorithm of the Revised Simplex 
Method, USAF Project Rand Memorandum, ASTIA Document No. AD 114135, 26 
Oct., 1953. 

eH M. Markowitz, The Elimination Form of the Inverse and Its 


Application to Linear Programming, Management Science, Vol. 3, No. 3, 
April, 1957. 
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The actual transformation may be summarized 


(1) Hay ] Km 1 2 .n 

(2 ) Lik 

(3) Lik EEST. 

(4) I: pour pP (j = r) kax41,2....n 


jk ” jk kr 

(5) lik = ER = u Lik (jr) 
which states (for those not versed in matrix notation), givenanxn 
matrix H (1), construct the corresponding n x n unit matrix I (2). 
Multiply the as column of H by the unit matrix (or the partly trans- 
formed matrix) to produce the pivotal column e (3). Divide the ele- 
ments of the pan row of the unit or matrix under transformation by Lih 
element of the pivotal column, (4). Correct all other rows by subtracting 
from the element to be corrected the product of pivotal element in that 
row and the corresponding element in the mE row computed in Step (4), 
(5). 


A simple example will make the procedure clear, say H is 


4 2 3 
102 
2 1 2| 


We start by writing down the unit matrix of corresponding order, viz., 


(2) 1 0 0 


0 1 0 
¡0 0 1 
K x DL. ROw T 
r 
(3) 1 0 0 4 Note: This is akin to merely augment- 
0 1 0 1 ing the unit matrix by H, for R= 1, 
0 0 1 2 
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K 21, R#l 
r 


: = = "9 9 6 
(09117570 0 4 Deal. c aatem 
0 12:0: 1 ] = r ぉ 1 
0 0 1 2 


TA AA OS AE 


Ih ¥ 2 m 


“1/4 1 O0 . Discard Pi . Note: Columns to the right 


j 


of i are not operated on, herein is the 
r 


great economy of this method. 


-1/2 0 1 


(3) 1/4 0 0 1/2 E . H, x P3 
“1/4 1 O -1/2 


-1/2 0 1 0 
(4) [1/4 0 0 1/21 Jo, * doy * E (second element in column P) 
E -2 0 -1/2 
-1/2 0 1 O | 
(5) 0 10. Iik = En z 2, D. 
MA ECL Note: Again since K = 2¢3, column 3 was 
-1/2 0 1 ・ unaffected; also note since P = 0, I 
2, 3k 
was unaffected. 
K = 3, our final transformation is: 
(3) 0 1 0 2 (^) [0 1 0 2 (5) Í 2 L -4 
-1/2 -2 0 -5/2 1/2 -2 0 -5/2 4/2 -2 5 
-1/2 0 1 1/2 -1 0 2 1/2 -1 0 2 


and therefore 


er 
a ey ee ln 
102 
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Check: 


[4] > 回 


[4 2 3 | | 2 l -4 1 0 0 
1 0 2 EE 5 - 0 1 0 
2 1 2 -1 0 2 0 0 l 


The total number of operations (neglecting all incidental nulls) is: 


Multiplications to bring in each new row n(n+1)(n+2)/3 
Multiplications and divisions to accomplish 
each transformation (n+l)(n)(n-1)/6 
To solve the equations once the inverse has 2 
been computed n 
Total product operations to complete solution nk + RE +n 
2 


Reference to Table III-l shows the efficiency of the elimination 
form of the inverse. It is important to note that by minimizing the 
number of operations, we also minimize the propagation of round-off error. 
Von Neumann has determined that an approximate inverse will be found 
if 


n «5 CN 


> where 3 is the base of the number system used in the 
computation, and S is the number of places, The program for the elimina- 


tion form of the inverse (appended) used a packed floating point where 


the word structure was: 


XS MMMM MMMM MQEE Legend: X = not used 
S = sign of magnitude 
M = magnitude (octal digit) 
Q = 2 binary bits of magnitude 
& sign of exponent 
E = magnitude of exponent 


=: Von Neumann and H. H. Goldstine, Numerical Inverting of Matrices 
of High Order, American Mathematical Society Bulletin, Vol. 53, pp. 1096. 


LIÉ 5 


E 





This realized an S of 29, and therefore theoretically the program should 
successfully compute the inverse for all matrices of n<23. 


The maximum round-off error to be expected is (after Von Neumann): 


-(S-1) /3 
と = .58 de 2a?) 


which for 3 = 2, S = 29, and n = 25 reduces to 


- 2 = .000067 


Thus we can expect an accuracy of, at worst, four significant figures for 
a 25 X 25 matrix. 

If we decide to use a double precision floating point, S (for the 
CRC 102A) would increase to 66, by using the following configuration: 

XS) MMMM MMMM MMMM XS, MMMM MMMM MMEE 


where Si is sign of magnitude 


S is sign of exponent 
and the round-off error now becomes (for n s 25) 


1.8125 -9 
€ * 3-7 x 1ol$ < 10 


This requires a considerably greater amount of computation time 
since a floating point multiplication is now three "old" floating point 
multiplications and two floating point additions, however even for n s 60, 
E EO For comparative purposes, note Table III-2. The choice of 
either packed or double precision is really dependent on the order of 
the matrix expected and the requírements of the problem, it seems to this 
writer that both procedures should be available. 


It must be mentioned in passing that the method of Hotelling and 


Bodevir- for iterating to find an improved inverse from an approximate 


1 
A. 5. Householder, Principles of Numerical Analysis, McGraw-Hill, 


New York, 1953, pp. 80 ff. 
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inverse and the original matrix, is readily available. However, although 
the programming of this scheme would be extremely simple, inasmuch as 

the basic matrix multiplication is already present, two A multiplica- 
tions are required for each successive improvement, and furthermore, 
another n x n matrix must be stored. These requirements with regard to 
storage and computation time make the double precision floating point 
arithmetic routine more attractive, since the time to solve an n x n and 
perform one iteration would be proportional to ont an + x 200 for 


the packed floating point, whereas with equal storage space, the double 


3 2 
precision would only be proportional to LA 425. 


The decision as to whether or not the simple packed floating point 
routine or the double precision was to be used would be made in the initial 
stages of the problem; i.e., when "n + r" was determined, assuming the 
number of elements was not specified (see Appendix I). 

In order to obtain a practical estimate of the accuracy of the 
matrix inversion program and the time required, the following tests were 
made: 

A symmetric 10 x 10 matrix was formed from a column of random numbers, 
the sign and position of the decimal point were also established using 
random numbers. This matrix was then inverted using the program. The 
matrix product of the original matrix and the inversion was then com- 
pared with a unit matrix. The rms. error of each element was found to 
be 431 x Kr en eco Se 

A similar test for 6 x 6 matrix resulted in an rms. error of 14 x 
9: 


10” and a maximum error of 43 x 10 


The running time was determined to be approximately .064 n minutes, 
(64 minutes for a 10 x 10, 14 minutes for a 6 x 6). 


ELE 37 





The flow chart and program for the elimination form of the inverse 


using simple packed floating point airthmetic is appended hereto. 
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TABLE III-1 


COMPARISON OF VARIOUS METHODS 


Method No. of Operations 
l 3 
Stiefel-Hestenes ZU does 
Strict ne I + Ae 
Peehovonalization. I + ae 
2 3 2 

Elimination Form 1/2 (n -2n +n) 
Back-Substitution 1/6 (ane + Te 


*Assumes product operation requires 70 word-times. 


TABLE III-2 


COMPUTATION TIME IN WORD-TIMES 


Operation Normal Float Packed Float 
Multiply 80 200 
Divide* 112 232 


*Assuming overflow occurs 1/4 of the time. 


(9) 
(24) 


ono del 


2 
Von Neumann, 
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Word-Times (n=#20)* 


112 x 10^ 


115 x 10^ 


113 x 10^ 


31 x 10^ 


2l x 10° 


Double Precision Float 
425 


684 





0123 
0124 
0125 
0126 
0127 
0130 
0131 
0132 
0133 
0134 
0135 
0136 
EF» 
0140 
0141 
0142 
0143 
0144 
0145 
0146 
0147 
0150 
0151 
0152 
0153 
0154 
0155 
0156 
0157 
0160 


MATRIX INVERSION PROGRAM 


INST 


m1 


0076 
0077 
0371 
3000 
3000 
0377 
0077 
0106 
0375 
2100 
0371 
0337 
0367 
0114 
0373 
0371 
0335 
0334 
3000 
0336 
2100 
2100 
0362 
0365 
2100 
0327 
3000 
0327 
2006 
0510 
2000 
3000 
1704 
3400 
0364 
3000 
2000 
0323 
0135 
0327 
0371 
0135 
0364 
0275 
0371 
3000 
0135 
0152 
3000 


M? 


2100 
2100 
0561 
3000 
2100 
2100 
2100 
0376 
0106 
2100 
0555 
0532 
2100 


0555 
0373 
2100 
0366 
2100 
0335 
0332 
2100 
0363 
0336 
2100 
0326 
2100 
0361 
0363 
2100 
2100 
2100 
2100 
2100 
2100 
2100 
2100 
0555 
0560 
0555 
0525 
0560 
2100 
0335 
0373 
2100 
0360 
0335 
2100 


M3 


0371 
0366 
0355 
0370 
2000 
0106 
D 
0106 
0106 
0373 
2006 
0114 
0514 
0114 
0373 
011) 
0336 
0365 
0124 
0336 
0152 
0373 
0135 
0327 
0323 
0133 
0145 
2006 
0140 
2000 


0145 
2001 
o400 
2001 
o400 
0364 
0323 
0135 
0327 
0131 
0135 
0000 
0375 
0156 
0161 
0135 
0152 
0127 


III ~- 1 - 1 


(CRC - 102A) 


REMARKS 


Store N 
Store (n +r) 
Shift N to M? position 
Clear cells 

0500 - 1777 
Initialize 0106 
Shift the 

contents of Channel O 

to Channel 17 
Set RC (row count) to 0 
N + 1 modifier 
Initialize 0114 
Build 

N x N 

unit 

matrix 
Set CM (colum count) to 1 
1701 + (n + r) = S 
Skip 0123 initially 
Increment CM by 1 
Set putaway to 0000 
Set RC to 0 
Reset pickup for ajy 
S - CM = RS( row selector) 
Clear CC (colum tally) and TS 
If RS greater 1677, skip to 133 
No, skip tp 0145 
Shift RS to M! position 
Intalize 0140 to RS 
Pickup a,, to "A" 
If A 4 O, skip to 0140 
If A= 0, skip to 0145 
load colum element 
Transfer to multiply 
load TS 
Transfer to add 
Store sum in TS 
Increment CC by 1 
Increment a 3 pickup by 1 
Increment aa by 1 
If N greater than CC, repeat 
No, compensate pickup 
Putaway row-column product 
Add 1 to RC 
If N greater than RC, skip 0155 
No, skip to 0155, thence to 0161 
Add 1 to 813 pickup 
Increment putaway address 
Repeat loop 





ADD 


0161 
0162 
0163 
0164 
0165 
0166 
0167 
0170 
0171 
0172 
0175 
0174 
0175 
0176 
0177 
0200 
0201 
0202 
0203 
0204 
0205 
0206 
0207 
0210 
021] 
0212 
0213 
0214 
0215 
0216 
0217 
0220 
0221 
0222 
0223 
0224 
0225 
0226 
0227 
0230 
0231 
0232 
0233 
0234 
0255 
0236 
0237 
0240 
0241 
0242 


MATRIX INVERSION PROGRAM 


INST 


36 
26 
35 
9 
20 


mM? 


0336 
0337 
2006 
2100 
2007 
0356 
0357 
2006 
2007 
0510 
0002 
3600 
2000 
0323 
0336 
3000 
0172 
0175 
3000 
2100 
0336 
3000 
0337 
0322 
0571 
05539 
0356 
2100 
0324 
0322 
0575 
03 36 
3000 
0510 
0002 
3,00 
2000 
0510 
3200 
2000 
0515 
0327 
0356 
0325 
3000 
0223 
0226 
0230 
3000 
0371 


M? 


0225 
0371 
0537 
2100 
0330 
0363 
03 30 
0363 
0332 
2100 
2100 
2100 
2100 
0335 
0323 
2100 
0360 
0335 
2100 
2100 
0325 
2100 
0325 
0330 
0335 
2100 
0363 
0363 
0363 
0332 
0336 
0375 
2100 
2100 
2100 
2100 
2100 
2100 
2100 
2100 
0335 
0375 
0360 
0335 
2100 
0360 
0355 
0371 
2100 
0336 


mM 


0357 
2006 
2007 
0323 
0356 
0172 
2006 
0175 
0175 
2000 
2001 
0400 
0510 
0323 
0201 
0204 
0172 
0175 
0172 
0325 
0207 
0242 
0522 
0324 
0327 
0373 
0222 
0223 
0226 
02 50 
0222 
0222 
0231 
2000 
2001 
0400 
2001 
2000 
OLOO 
0510 
0373 
0236 
0356 
0325 
0205 
0223 
0226 
0230 
0217 
0123 


111 - 1 - 2 


(CRC = 102A) 


REMARKS 


CM = ] = J 
J xN 
0500 + (J x N) 
Clear CC 
Initialize Iri 
MNO (ご 
Initialize Py 
in 0173 E 
Initialize It 
Load T 
Load Py 
Ir x Py = i 
Putaway result 
Increment CC by 1 
If CM greater than CC, 0201 
If CM = CC, 0204 
Increment I to I 
Increment nn, ried 
Continue dividing row 'r' 
Set M to O 
1f CM greater than M, 0207 
If CM = M, 02X2 
Build I', 
Build Ty 3 
N + 1 to TS 
Initialize RC to 1 
Initialize Iip pickup 
Initialize Py pickup 
Initialize I’! pickup 
Initialize I; , pickup 
If RC 4 CM 
skip to O222 
If RC = CM, skip to 0231 
Load oe 
Load B, 
Multip1F 
Shift for subtraction 
Load 1 


in 0175 


i. = =) jak or 
Put bay E F Er 
Increment RC by 1 

If N + 1 greater RC, 0236 
Increment I), to Ia, 
Increment M by 1 

Repeat loop 

Increment Per to Pre] 
ox + N = I5 

Increment nut 

Repeat loop 

Next iteration ( « loop) 





ADD 


0243 
ooh 
0245 
0246 
0247 
0250 
0251 
0252 
0253 
0254 
0255 
0256 
0257 
0260 
0261 
0262 
0263 
0264 
0265 
0266 
0267 
0270 
0271 
0272 
0273 
0274 
0275 
0276 
0277 
0300 


0301 


0322 
0323 
0324 
0325 
0326 
0327 
0330 
0331 
0332 
0333 
0334 
0335 
0336 
0337 


MATRIX INVERSION PROGRAM 


INST 


25 
>> 
30 
32 
>> 
3h 
3D 
うら 
Do 
52 
3h 
の の 
う フ 
32 
> 
30 
32 
34 
>> 
36 
30 
30 
3h 
30 
34 
30 
34 
30 
34 
52 


- 0521 


M+ 


2100 
0347 
0367 
0376 
0354 
3000 
0346 
0353 
2100 
0351 
3000 
0352 
0276 
0351 
a_100 
0352 
2006 
3000 
0267 
0270 
1701 
0000 
3400 
0364 
3000 
2000 
0270 
0364 
0366 
2100 


M? 


2100 
2100 
2100 
0332 
2100 
2100 
2100 
2100 
2100 
0332 
2100 
0435 
0435 
0363 
2100 
0361 
0363 
2100 
0360 
0360 
2100 
2100 
2100 
2100 
2100 
2100 
0350 
2100 
0352 
2100 


M? 


0556 
0152 
0000 
0152 
0155 
0125 
0152 
0195 
0352 
0276 
0260 
0352 
0276 
0267 
0364 
2006 
0270 
0267 
0267 
0270 
2000 
2001 
OLOO 
2001 
OLOO 
0364 
0265 
1001 
0256 
0345 


(CRC - 102A) 


REMARKS 


Set CM = 0 

Reset putaway 

Putaway "1" in 0000 

Dummy (see 02414) 

Mod 0155 ("out com") 
Transfer for N+1°° colum 
Reset putaway 

Restore 0155 

Set EC (equation count) = 0 
Initialize putaway 

Skip modifiers 

Increment EC by 1 

Increment putaway 

Reset ho 

Clear TS 

Shift EC 

Extract into b pickup 

Jump modifiers 

Increment h 

Decrement b 

Load h 

Ioad b 

Multiply 

Load TS 

Add 

Putaway summation 

If b greater than O, repeat 
No, summation to answer storage 
If (n + r) greater EC, 0256 
Clear tally, inversion completed 


Available for boot-strapping routines 


Denominator in Channel O 
Numerator in Channel 1 


0000 
0000 
0502 
0000 
0000 
0000 


0000 
0000 
0000 
0000 
0000 
0000 
0000 
2005 
0000 
2100 
0000 
0000 
0000 
0000 


0502 
0003 
0000 
0003 
1677 
1705 
0030 
2000 
IF 
0105 
1701 
0001 
0000 


0500 


III-1-3 


I';, in M storage 

CC storage 

D M! storage 

M ¿torage 

Testword for 0131 

RS storage 

Shift Control 

Clear Routine (2001) 
extractor 

Clear Routine (2003) 

Constant for 0116 

Constant for Clear Routine 

CM storage 

Constant for 0113 etc. 





MATRIX INVERSION PROGRAM 


ADD INST 
0340 00 
0341 00 
0342 10 
0343 00 
O344 00 
0345 00 
0346 30 
05,7 36 
0350 00 
0351 00 
0352 00 
0555 34 
0354 34 
0555 00 
0556 00 
0557 00 
0360 00 
0361 OO 
036z 00 
0565 00 
0364 02 
0365 00 
0366 00 
0367 00 
0370 26 
0571 00 
0372 3, 
0575 OO 
0374 00 
0375 00 
0576 00 
0577 30 
Channel kl 
0400 30 
0401 32 
0402 30 
0403 32 
ohoh 27 
OhO5 32 
0406 32 
0407 31 
OL10 32 
0h11 27 
0h12 52 
0h15 32 
Ohih 36 


Mi 


0000 
1000 
0100 
0000 
0000 
0000 
0364 
2100 
0000 
1700 
0000 
3000 
3000 
0003 
0511 
0000 
0001 
0000 
0500 
7777 
2314 
0000 
0000 
4000 
2100 
0000 
2004 
0000 
2100 
0077 
0001 
0000 


3000 
2006 
2001 
2000 
2002 
2000 
2000 
2002 
2007 
2002 
2007 
2007 
2006 


M? 


0000 
2100 
2100 
0000 
2100 


contains packeä 


0443 
0475 
2100 
043% 
0451 
0453 
Ok 44 
2003 
Ol 5l 
0451 
0436 
Obl 
0435 


M? 


0000 


floating point 


2006 
0433 
2007 
2102 
2002 
2002 
2103 
2000 
2102 
2002 
2002 
2105 
2006 


III-1-h 


(CRC = 102A) 


REMARKS 


Available 
Constant for 0156 
Not required 

Not required 

Not required 

Not required 
Reset for 0152 
Reset for 0244 
Testword for 0275 
Putaway reset 

EC tally 

Modifier for 0252 
Modifier for 0247 
N in M? 

I in M? 

J tally 

Constant for 0146 
Shift control word 
Constant for 0126 

M! extractor 
Temporary storage (TS) 
S storage 

(n + r) storage 

"1" in PackBinFitPt 
Clear Routine (2000) 

N storage 

Clear Routine (2002) 
RC tally 

Clear Routine (200) ) 
Testword for 0105 
Constant for 0104 

Not required 


arithmetic routines 


Prepare 

exit 
Store 204 operand 
Extract sign of exp? 
Shift to operating location 
Extract mag of exp! 
Extract sign and mag of mag 
Shift to 2000,2001 
Extract sign of exp? 
Shift to operating location 
Extract mag of exp? 
Extract sign and mag of mag? 
Pickup entry command 


1 





ADD 


0415 
0416 
0417 
o420 
0h21 
0122 
0423 
ohol 
0425 
0426 
0427 
0430 
0431 
0432 
0433 
0434 
0435 
0436 
0437 
OLLO 
0441 
Oke 
0443 
0444 
OLS 
0446 
0447 
0450 
0451 
0452 
0453 
0454 
0455 
0456 
0457 
0460 
0461 
0462 
0463 
046% 
0465 
0466 
0467 
0470 
0471 
0472 
0473 
Oh Th 
0475 
0476 
0477 


MATRIX INVERSION PROGRAM 


INST 


30 
32 
32 
34 
Do 
36 
30 
35 
37 
3] 
27 


M1 


2006 
2006 
0273 
2007 
2002 
2002 
2003 
2001 
2001 
2000 
2002 
2002 
2004 
2003 
3000 
0000 
0000 
0000 
0000 
2007 
2100 
3000 
0000 
Wart 
2007 
2000 
2001 
3200 
0000 
0000 
0000 
2000 
2001 
2001 
3000 
0000 
3777 
2001 
2001 
2001 
2000 
3000 
2000 
2001 
2002 
2003 
2001 
3400 
0000 


0000 
0000 


M? 


0437 
0461 
0h61 
O&OO 
2000 
2000 
2006 
2007 
3000 
2001 
0h52 
0436 
0476 
2100 
2100 
0000 
0000 
0000 
0000 
0h50 
2005 
2100 
0000 
ferma 
047% 
2002 
2003 
2100 
0000 
0000 
0000 
2002 
2003 
3000 
2100 
0000 
0000 
0460 
0460 
0477 
0477 
2100 
2002 
2006 
2100 
2007 
3000 
2100 
0000 


0000 
0000 


M? 


2006 
0h17 
2107 
OL4O 
0467 
2006 
2007 
2001 
0462 
2002 
2104 
2004 
2003 
2000 
0274 
0100 
0001 
0077 
0030 
0445 
2003 
0h21 
0030 
7600 
0454 
2000 
2001 
0426 
0037 
0037 
0077 
2000 
2001 
0462 
0526 
0001 
0000 
2001 
2001 
2001 
2000 
0426 
2006 
2007 
2000 
2001 
0,62 
0426 
TITI 


0177 
0001 


PIT OS う 


(CRC = 102A) 


REMARKS 


Shift to M! 
Set 0417 
Extract control 
If control > 3000, O&O 
No, operations 
for normal 
floating point 
addition 
routine 
Store in 2002, 2003 
Shift sign of exp for pack 
Extract mag of exp for pack 
Pack exponent 
Packed result to 2000 
Exit to main program 
Extractor for sign of exp 
Constant for 0414 
Extractor for mag of exp 
Shift control word 
If control > 3200, 0445 
Change sign of 2P% operand 
Transfer to addition 
Shift control word 
Extractor for magnitude 
If control a 3400, 0454 
No, operation for normal 
floating point 
multiplication 
Shift to pack commands 
Shifter for sign of exp 
Extractor for mag of exp 
Operations for 
normal floating point 
division 
Transfer to pack commande 
Constant for 0462 
M! extracotor 
Normal 
correct 
overflow 
procedure 
Transfer to pack commands 
Operations for 
normal floating 
point addition 
when 2002 > 2000 
see O42) 
Transfer to pack commands 
M? extractor 


Sign and exponent extractor 
Constant 











Clear 0500-1777 
[Hj— 1700 

set RC = 0 
Set N+l in modifier 






This N is the 
order of the 
matrix, viz. n 


Initialize "1" PA 
















Add N+l to 
u PA 


Command 


a 


13612) 
add l to RC 








Unit Matrix 





Set 1—0CM 
1701 +n + r>S Be 
© "Add 1 to CM - 
[CM] -1=> J 
| Add 1 to a, Jj*N + 0500 = n 
Reset PA to 0000 | Pickup 0000 « [j) - P, 
up to ¡Add 1 to putaway z 
Set RC to 0 P P 
multiply Set CC to O 
Reset a to a in rows 


1 pickup 


1j l 





S - CM>RS 
Clear CC 
Clear TS 


Row selection 















Yes 





A]. RSH B 
[ts] + [J>TS 






Add 1 to CC 


No 


RS + 1—> RS 


Add 1 to a 


1 1 pickup 





Putaway 2 TS — 0000 
(initially) 
Add l to RC 
























d 
1, + P? Ii 


r r r 
Add l to CC 











Mod to l, + M I, 


r r r 






















Add l to CM => CM 


Reset "ij pickup to a 


PA "1" NPFP—> 0000 


Set PA to 0001 
Alter "out" cmd to rtn to ® 







11 











Yes 













0500 + M—> [I 
j 


Set P—> 0000 
Y 


Set Dec to N41 


Set RC to 1 











A 
Restore "out" cmd 
Set EC to 0 
Yes No set putaway 
Reset 
CMDR add 1 to EC "EC 
Clear TS 






Mod PA Shift EC 


Initialize b. pickup 







TS + (hic, )* (b,) — TS 






No 





b,» TW 






IW = 000021002001 






put away 


$ Kote © 








APPENDIX IV 
LAGRANGIAN INTERPOLATION METHOD FOR 


THE SOLUTION OF ALGEBRAIC EQUATIONS! 


The Lagrangian interpolation method described by D. E. Muller was 
selected for the solution of large degree polynominals. This procedure 
was programmed, with certain modifications, and tested on the CRC-102A 
computer. By removing the roots, as found, by synthetic division, rather 
than by the algorithm recommended by Muller, accuracies of at least eight 
decimal places were obtained for all roots of equations up to degree 21. 
This was considerably better than the accuracies reported by Muller, and 
incidentally, probably entailed considerably more iteration. 

A graph of solution time vs. degree of equation, for € = ie is 
presented in Figure IV-1. Convergence to a root generally requires be- 


tween 7 and 14 iterations (for the € of lo The iterative loop is 


completed in about 1.5 minutes, varying slightly (not more than 107) 
according to the degree of the equation. In making these tests, it was 
found that when equations had roots of the same modulus, there is a pro- 
nounced slow down in the rate of convergence, The equations xd +1=0, 
and Roy + 1 = 0 required about 12.5 minutes per root, whereas the more 
general polynominal, (equation A under Results of Tests) only required 
about 9.4 minutes per root. Round-off error was not present in the 

20 24 


x + 1 equation and affected only two roots in the solution of x +12 0 


within the eighth decimal place. 


=D. E. Muller, Solving Algebraic Equation by an Automatic Computer, 


Mathematical Tables and Other Aids to Computation, Vol. X, No. 56, Oct., 
1956, pp. 208-215. 


je 1l 





Thus far the program has failed only once. When solving an equation 
with a real double negative root and a pair of complex roots of about 
the same modulus, viz., 4; the rate of convergence was normal until about 
four decimal places had been established and decreased so as to become 
barely perceptible. The program did, however, remove all roots interior 
to modulo 4 without difficulty. 

The basic outline of the Lagrangian interpolation method is as 
follows: 

Successive iterations toward any particular root are obtained by 
finding the nearer root of a quadratic whose curve passes through the 
points; ZU £(2,); Zi]? f(z,_,); Zio? £(z, 5). 

The Lagrangian formula may be written as L. [ gs) | = bo + biz * b, 
where the b's satisfy the system of equations: 

2 


by? + bizi + b, = £(z,) 


De + bu wel, 7) 


=] 


by? 


o 
eN eN 


> + b121.» a b, = f(z,_,) 

and these equations, of course, become almost one in the same as the 

iterative process converges, depending on the accuracy which we specify. 
We are to construct a polynominal (quadratic) whose curve passes 

through the three points. By the Lagrangian interpolation ee we 


write: 
(zw mz ze) (ZZ zu.) 
EPA) AA M 
E ES NU Cea ce pr? Gp up XML 


CES TCU D 





a CS EE D 


zn E. Milne, Numerical Calculus, Princeton University Press, Princeton, 
N. J., 1949, for example. 


IV = 2 


x 





If we make the following simplifying substitutions: 
h/h, = A 


El =) 


L+A = A. 


we may express the coefficient of f(z,) as: 
2 -1 -1 
NA gt AAPOAD +1 
Dag uu. 
7) 


UA. 


(CREE e.) 


since 
AT 2) 


i i 
1 Ber er N 
l + x + * = z gem = 5 + N 
T i i-l i 
BA: Zz = Z 
A bv | = TIA and finally 
Aj ee 
VIN [Ai A E 
AN N > y See «ss; therefore 
i} ry gue 513 | 


z. jJ(z 7 2z, 5) 
= i-2 the coefficient of f(z); 


(z - 
) 


-1 
e) A MAD ) = (@ = 2 )@,- 2,5 


which when expanded is: 


\ Ai A thr, Aj + 入 + ュ = 人 N N e AA, +A) + 1 (QED). 


By similar manipulation, the coefficients of f(z,_,) and f(z, 7) can be 


found and the desired quadratic in ut 


IV - 3 





ET ES) - Aj Om A? peo a4 AA, S D | + 


-1 2 2 ] 
AA [ce f(z,_,) an + NUS +A) + f(z,) 
A single iterative step is obtained by letting 2 be such a value 


of z that L; [f(z] z 0. We may solve the quadratic inversely to obtain 


= Z 2 


Z 
which is A = tht from which we obtain the new z, and repeat 
i+l Z =Z] i 


the iterative process. 


All operations are in complex numbers and in floating point arithmetic. 


A root is taken when 


a; 
Z 


i < P 








For proof of convergence of the process the reader is referred to Dr. 


Muller's paper. 


IV - 4 


RESULTS OF TESTS 


Equation A 
ee ee 264207215 .- 479,0 2t 
EU E c15405527 79782002 20890122500 z^ + 912.75 2° 


+ 254.75 = + 1411.0 z^ + 1319.0 z^ + 1022.5 の + 470.5 z 


Pos De LO >) 


Total solution time: 2 hours, 34 minutes, 10 seconds on CRC-102A. 


Roots: 
-0.284460054 + j 0.598759363 
-0.376766336 + j 0.188840569 
-0.101279181 + j 1.429775042 
-0.0449 75589 + j 1. 738162299 
+1.999999996 
-0.011258189 + j 1.933289317 
+2.500000009 
-0.181260648 + j 1.034830851 
+3.000000001 
*3.499999991 


IDEs 


Equation B 


30 


RESULTS OF TEST (Continued) 


x + l = 0 


Total solution time: 


Roots: 


-0.000000000 
-0.951056516 
+0.866025403 
-0.587785252 
+0.994521895 
-0.207911690 
-0.994521895 
+0.40673663 7 
+0.207911688 
+0.743144829 
-0.743144826 
+0.587785258 
-0.86602 5404 
+0.951056517 


-0.406736643 


E = 2 


6 hours, 23 minutes, 30 seconds on CRC-102A 


+ 7 I 1 IF + d* d r dg è は は p 
に 


{+ 
C. 


IV - 6 


Dr 


0.999999999 
0.309016994 
0.500000000 
0.809016994 
0.104528463 
0.97814 7601 
0.104528463 
0.913545463 
0.978147598 
0.669130603 
0.669130607 
0.80901 7000 
0.499999998 
0.309016993 


0.913545458 


13 


) 


6 








の 
= 
MQ 
Eh 
® 
> 
qo 
= 
ord 
E 
8 
o 
ord 
42 
2 
rd 
O 
u) 

















Degree of Equation 





pur 
FF 








+ oon EBEG 
E tB 
i キー トキ E. LES 


a hot + ELA ー 
des 


ー ォ ーー ペー オー A W ュー っ = 一 


u uot nio 




















LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD 


INST 


M+ 


M? 


MS 


REMARKS 


Channels O aná 1 are used temporary storage. 


Channel 2 contains a decimal to binary floating point conversion 
which while not needed in overall program is included here. 


0200 
0201 
0202 
0203 
0204 
0205 
0206 
0207 
0210 
0211 
0212 
0213 
o2ik 
0215 
0216 
0217 
0220 
0221 
0222 
0223 
0224 
0225 
0226 
0227 
0230 
0231 
0232 
0233 
0234 
0235 
0236 
0237 
0240 
0241 
0242 
0243 
0244 
0245 
0246 
0247 
0250 
0251 
0252 


26 
3D 
52 
55 
52 
05 
30 
32 
30 
30 
3) 
xl 
32 
52 
26 
うう 
26 
30 
3) 
31 
59 
34 
ah 
00 
00 
5h 
34 
00 
z2 
32 
25 
30 
52 


1340 
0246 
0251 
0255 
2100 
3000 
0254 
0254 
0016 
0017 
0230 
3000 
0230 
0230 
2004 
2006 
2004 
0230 
0230 
2000 
2000 
0247 
3000 
0000 
0000 
0247 
3000 
0016 
0247 
0247 
2003 
0265 
0247 
2003 
0260 
2001 
2002 
5000 
0002 
0000 
0002 
0000 
0000 


0250 
2005 
2100 
2100 
0252 
3000 
2100 
0255 
2100 
2100 
2100 
2100 
0255 
0265 
2005 
2001 
0257 
0262 
2100 
2001 
0261 
2100 
2100 
0000 
0000 
2100 
2100 
2100 
0255 
0265 
0263 
0256 
0265 
0263 
0265 
2100 
2003 
2100 
2100 
0000 
0000 
2100 
0000 


2005 
0255 
0210 
0211 
0270 
5000 
0265 
2004 
0230 
0247 
0214 
0231 
0227 
2005 
2006 
2001 
2004 
0230 
0215 
2000 
2000 
0235 
0267 
0000 
0000 
02 34 
0267 
0230 
0227 
2003 
2003 
0265 
2003 
2003 
0237 
0266 
2000 
0267 
0230 
0000 
0000 
0250 
TTTT 


IV - 1 - 1 


Build and 
store testword 
Reset pickup 
Reset pickup 
Reset putaway 
Clear buffer 
Reset extractor 
Set 2004 to 1 
Pickup integer 
Pickup fraction 
Is integer > 0? 
No, test fraction 
Yes, extract sign 
Extract least sig. digit 
L.S.D. x (2004) 
Accumulate products 
Increase 2004 by power of 10 
Shift decimal word hà right 
If conversion not complete, 0215 
If complete normalize 
and correct exponent 
Is fraction > 0? 
No, skip to putaway 
Sign storage 
Integer storage 
I = O, Is fraction > 0? 
If both are O, putaway O 
Testword 
Extract sign when fraction only 
Extract least sig. frac. dec. dig. 
Divide by 10 
Shift extractor 
Extract next digit 
Divide by 10 
If conversion not complete, 0257 
If integer àa O, transfer to add 
If integer was O, normalize 
and putaway 
Testword addend 
Fraction storage 
Testword build constant 
Integer pickup reset 
M* extractor 











LAGRANGIAN INTERPOLATION METHOD FOR HIGR DEGREE POLYNOMIALS 


ADD 


025% 
0254 
0255 
0256 
0257 
0260 
0261 
0262 
0263 
0264 
0265 
0266 
0267 
0270 
0271 
0272 
0273 
0274 
0275 
0276 
0277 


INST 


M1 


0001 
0000 
0000 
0000 
0000 
7400 
0000 
0000 
5000 
0000 
7400 
3000 
0227 
2000 
c270 
2007 
2001 
0210 
0211 
0270 
0233 


M? 


2100 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
2100 
0255 
2100 
0255 
0252 
2100 
0250 
0250 
0264 
0210 


mM 


0247 
0017 
0001 
000% 
0012 
0000 
0044 
0004 
0000 
0002 
0000 
1540 
2001 
0016 
2007 
0273 
0015 
0210 
0211 
0270 
0205 


REMARKS 


Fraction pickup reset 
Extractor reset 
Sign extractor and -1 
Shift control - h left 
"Power of 10" 
Extractor testword 
To correct exponent 
Shift control = 4 right 
10 for division 
2 
Extractor working storage 
Transfer to addition 
Extract sign into mag portion 
Putaway exponent 
Build magnitude 

putaway 
Putaway magnitude 
Modify integer pickup 
Modify fraction pickup 
Modify putaway 
Have all numbers been converted? 


Channels 3 through 7, and 1000 through 1012 contain main program 
for Lagrangian Interpolation Method; in the program detailed 
here Channels 10 and 11 contain & binary floating point to 
decimal conversion routine which would not be required for 


the overall program. 


Channel 15 is used for tempoary storage 


and Channels 14 through 16 contain the floating point arith- 


metic routines for complex numbers. 


IV -1-e2 











LAGRANGIAN INTERPOLATION METHOD FOF HIGH DEGREE POLYNOMIALS 





ADD INST M1 M? M3 REMARKS 

0300 30 1340 2100 1372 N to N* 

0301 05 3000 3000 0000 Transfer 

0302 Ok 3000 3000 0100 the contents 

0205 735 050] 1375 0501 of Channel O 

0304 45 0302 1376 0302 to 

0305 3h 1374 0302 0301 Channel 1 

0306 36 0301 1437 0301 and 

0307 36 0302 1437 0302 reset, 

0310 30 re. 2100 2000 

050T 36 2100 1907 2001 Initial setup 

0312 31 2000 2001 1500 for 2; _2 

2515. 52 155 Tf 1520 

0314 321 1331 Bist 1301 Initial setup 

05315 31 1533 1337 1321 for 24.1 

05316 31 1353 1337 1302 Initial setup 

0317 351 1333 1337 1322 for 2; 

0420 45 2100 2100 1303 Initial setup 

0321 36 2100 1337 1314 for Ai 

0522 31 1333 1337 1323 - 1/2 + jo 

0523 26 1372 1466 2007 Commence pickup build 

0324 26 2100 2100 1434 Clear 1434, 1475 

0325 26 2100 2100 1435 Clear 1435, 1476 

0326 35 0640 2007 0327 Set pickup for A 

0327 31 0102 0103 2000 Pickup Ao 0 

0530 31 2000 2001 1507 Ao 1s (eof f(zi) 

$551 51 135 037 1327 daf (zi) = O 

0332 36 0527 1473 2007 Build Ag pickup 

0555 42 2007 1464 0334 and enter in 0354 

0334 41 0070 0071 2002 A, pickup command 

0535. 34 3000 2100 1540 Transfer to floating addition 
0336 36 0335 1473 0334 Modify A, pickup for Al. 
0337 3h 0640 0334 0341 If sum complete, putaway 
O340 34 3000 2100 0334 No, repeat to 0334 

0341 34 1475 2100 0347 Has even sum putaway been made? 
0342 31 2000 2001 14 34 No, even sum to 1434, 1475 
0543 36 0327 166  0O3kl Build Aj, pickup 

OSshh 31 0100 0101 2000 A pickup 

0345 36 Oshh 14735 2007 Set An pickup for odd terms 
OF46 43h 3000 2100 0555 and repeat back for odd sum 
0347 32 2000 2001 1425 Putaway odd sum in 1435, 1476 
0350 31 14 34 1475 2002 Load even sum 

0551 34 3000 2100 1540 Odd + even for 

0352 31 2000 2001 1306 f (24.1) 

0355 31 f 1537 1526 and clear imaginary 

03584 31 1434 1475 2000 Load even sum 

0355 31 1435 1476 2002 and odd sum 

0356 3% 3000 2100 1516 for subtraction 

0357 31 2000 2001 1305 to form f(z;.9) 








LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD 


0360 
0361 
0362 
0363 
0364 
0365 
0366 
0367 
0370 
0371 
0372 
0373 
0374 
0375 
0376 
0377 
OLOO 
OL01 
0402 
0403 
OL OL 
0405 
0406 
0407 
0410 
0411 
0412 
0413 
0414 
0415 
0416 
0417 
0420 
0421 
0422 
0423 
oh ol 
0425 
0426 
0427 
0430 
0431 
0432 
0433 
043% 
0435 
0436 
0437 


INST 


M2 


1533 
1303 
1323 
1331 
1533 
3000 
2000 
2002 
1306 
1326 
3000 
3000 
1305 
1325 
1303 
1323 
3000 
3000 
1307 
1327 
3000 
1550 
1332 
3000 
1304 
1324 
3000 
1303 
1523 
3000 
1307 
1327 
2004 
2006 
3000 
3000 
3000 
1304 
1324 
3000 
3000 
3000 
1303 
1323 
3000 
3000 
1303 
1323 


M? 


[557 
1344 
1364 
1555 
21551 
2100 
2001 
2003 
1347 
1367 
2100 
3000 
1346 
1366 
1344 
1364 
2100 
3000 
1350 
1570 
2100 
15372 
1373 
2100 
1345 
1365 
2100 
1344 
1364 
2100 
1350 
1370 
0264 
0264 
2100 
3000 


1345 
1365 
2100 
3000 
3000 
13,4 
136) 
2100 
3000 
1344 
1364 


M3 


1325 
2000 
2002 
2004 
2006 
1410 
1304 
1324 
2004 
2006 
1607 
1330 
2000 
2002 
2004 
2006 
1607 
1310 
2004 
2006 
1410 
2004 
2006 
1416 
2004 
2006 
1607 
2004 
2006 
1607 
2004 
2006 
2004 
2006 
1607 
1440 
1330 
2004 
2006 
1607 
1330 
1310 
2004 
2006 
1607 
1310 
2000 
2002 


V-1-4 


REMARKS 


end clear imaginary 
Load 

A1 
Load 

1 + jo 
Transfer to cpx. add. 
Putaway 

By 
Load 

Dzs-1)- 
Transfer to cpx. mpy. 
Store in B 
Load 


f(zı.2) 
Load 
Ai 


Transfer to cpx. mpy. 
Store in A 
Load 
f(z,) 
f(z E + f(z,) 
1 i 
l2 = の 


wiz; 
Fear, | -f(najA, * fU 


cz. ‘Ay 
Eu 


gil. ‘Ay A, 
Load 

f(z, ) 

ana, "multiply 


br(2,) 0 A, £(Z) = WeAC 


1 
Store in C 


Load 


A 
ia 
Store in B 
f(z1.2) Ay fromA 
Load X; 


Multiply for f(z}_2) d4* 
Store in A 
Load 


N 





LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD 


O40 
0441 
04402 
0443 
0444 
0445 
0446 
0447 
0450 
0451 
0452 
0453 
0454 
0455 
0456 
0457 
0460 
0461 
0462 
0463 
046% 
0465 
0466 
0467 
0470 
0471 
0472 
0473 
0474 
0475 
0476 
0477 
0500 
0501 
0502 
0503 
0504 
0505 
0506 
0507 
0510 
0511 
0512 
0513 
0514 
0515 
0516 
0517 
0520 


INST 


M+ 


1304 
1323 
3000 
1307 
1327 
3000 
1310 
1312 
3000 
1330 
1332 
3000 
3000 
2000 
2002 
3000 
1h50 
1442 
3000 
3000 
2000 
2002 
2100 
3000 
3000 
2000 
3000 
1434 
2000 
1440 
3000 
2000 
1331 
3000 
2000 
3000 
24354 
3000 
2000 
1435 
1331 
3000 
2000 
3000 
1434 
3000 
2000 
2000 
2000 


M? 


1345 
1365 
2100 
1350 
1370 
2100 
1351 
15545 
2100 
2571 
1319 
2100 
3000 
2001 
2003 
2100 
1401 
1403 
2100 
3000 
2001 
2100 
2003 
2100 
2100 
2001 
2100 
1475 
2001 
1401 
2100 
2001 
1337 
2100 
133] 
2100 
1475 
2100 
2001 
1476 
1337 
2100 
1331 
2100 
1475 
2100 
2001 
2100 
2100 


M? 


2004 
2006 
1410 
2004 
2006 
1607 
2004 
2006 
1410 
200% 
2006 
1416 
1330 
2004 
2006 
1607 
2004 
2006 
1416 
1440 
2004 
2006 
2007 
1607 
1551 
1434 
1551 
2004 
1434 
2000 
1507 
1435 
2002 
1540 
2000 
1551 
2004 
1500 
1436 
2002 
2000 
1516 
2000 
1551 
2004 
1500 
1434 
1451 
1453 


IV ~ 1 = 5 


REMARKS 


Load 
A 


i 
(Ay +4,) 
Load 


a (a A.) 
f ( ( 
a 

DAS 


er 


ai = E EAS * f(z, )(Ai+0) 
Result = 8,5; Store in B 
Load B for complex 
multiplication, i.e. B= 
Load AC 
Be = WAC 
Store in C 
Load 
BS = VAC 
for conjugate 
multuplication 
(B2 -hAC)(B? - hAC)" 
Compute r 
Store r in TS, 
Compute ir 
Load r 
Store Jr in 18, 
(B? = LAC) = 
x/r = cos 6 
Store cos 6 in TS, 
Load 1 
l + cos 8 
(1 + cos 6 )/2 
{MA + cos Q 7 = cos 9/2 


Load ir 


(1 - cosO)/? 
l =- cos O)/?2 = sin 6/2 
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ADD 


0521 
0522 
0525 
0521 
0525 
0526 
0527 
0530 
0531 
0532 
0555 
0534 
0555 
0536 
0557 
0540 
0541 
0542 
0543 
0544 
0545 
0546 
0547 
0550 
0551 
0552 
0553 
0554 
0955 
0556 
| CONT 
0560 
0561 
0562 
0565 
0564 
0565 
0566 
0567 
0570 
0571 
0572 
0573 
0574 
0575 
0576 
0577 


M1 


1436 
1436 
1475 
2007 
2100 
1477 
1403 
2007 
2100 
3000 
2007 
2100 
3000 
1450 
1451 
3000 
2000 
2002 
2100 
5000 
2000 
3000 
1452 
1555 
3000 
2000 
2002 
2100 
3000 
2000 
3000 
14 3) 
1436 
1475 
1451 
1450 
3000 
1453 
1452 
3000 
3000 
1307 
1327 
1304 
2100 
1324 
2100 


M? 


2100 
2100 
1575 
2100 
2007 
1515 
2100 
2100 
2007 
2100 
2100 
2007 
3000 
1460 
1461 
2100 


MS 


1450 
1452 
2107 
1461 
1463 
2107 
0790 
1462 
1460 
0535 
1460 
1462 
1330 
2004 
2006 
1416 
2004 
2006 
2007 
1607 
1434 
1330 
2004 
2006 
1416 
2004 
2006 
2007 
1607 
1436 
1330 
0566 
0563 
0566 
2006 
2004 
0570 
2006 
2004 
1410 
1330 
2000 
2002 
2004 
2005 
2006 
2007 


TV 1«6 


REMARKS 


PutawBy X 
and X, exps 
Set ¿o mag positive 
and 
A, mag negative 
Clear mag sign bit 
If Q.(B? - hAC) » O, 0553 
No, putaway “y positive 
and x, negative 
Skip to 0535 
Putaway «o9 positive 
and *, negative 
Load B4 
Load xo 
Load (Go 
Subtract 
Load for 
conjugate 
multiplication 
Multiply and 
store in TS, 
Load 8 
load x, 
Load (31 
Subtract 
Load for 
conjugate 
multiplication 
Multiply and 
store in TS, 
Load 84 
If TS, > TS,, 0566 
1f TSg > TS,, 0563 
If TS, > TS}, 0566 
Load x, 


load e 
and (2 
B + pe “LAC 
Store in B 
Load 
f(z, ) 
Load 
SEN 
for 1 
-2 Ar f(z, ) 
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ADD 


0600 
0601 
0602 
0603 
0604 
0605 
0606 
0607 
0610 
0611 
0612 
0613 
0614 
0615 
0616 
0617 
0620 
0621 
0622 
0623 
0624 
0625 
0626 
0627 
0630 
0631 
0632 
0633 
06 34 
0635 
0656 
0657 
0640 
0641 
0642 
064 3 
O6hl 
0645 
0646 
0647 
0650 
0651 
0652 
0653 
0654 
0655 
0656 


0660 


M1 


5000 
1520 
1332 
3000 
3000 
1302 


3000 
1330 
1332 
3000 
1302 
1422 
3000 
3000 
1302 
1522 
3000 
1302 
1322 
3000 
1454 
3000 
145% 
1372 
1470 
1470 
2007 
1310 
1312 
0100 
1555 
0657 
3000 
2002 
0677 
0104 
3000 
1602 
1310 
1312 
0646 
0646 
3000 
0104 
0000 
3000 


M? 


2100 
1371 
1975 
2100 
3000 
1343 
1363 
1342 
1362 
2100 
1371 
1373 
2100 
1343 
1363 
2100 
3000 
1343 
1363 
2100 
1343 
1363 
2100 
2000 
2100 
2002 
1466 
2007 
2100 


1351 
1555 
0101 
1557 
1636 
2100 
2003 
1636 
0105 
2100 
1643 
1351 
13353 
1466 
0656 
2100 
0105 
0000 
3000 


M3 


1607 
2004 
2006 
1650 
1330 
2000 
2002 
2004 
2006 
1416 
2004 
2006 
1607 
2004 
2006 
1410 
1310 
2004 
2006 
1416 
2004 
2006 
1650 
0631 
0632 
0725 
2007 
0656 
2007 
0646 
2004 
2006 
2000 
2002 
1635 
1611 
1602 
1550 
2002 
1542 
2002 
2004 
2006 
0646 
0660 
1611 
2000 
0644 
1,40 


Dee 


REMARKS 


Multiply 
load 

B + „Be -LAC 
Divide for \ 
Store in B AO 
Load 


Load 


eet 
Zi = Z1-] - h, 
Load 


Ay 
Ay (ry) ln 


21 
24 * hi4) 7 244) 
Store in A 
Load 
24 
“ity io Ji 
ht EAP O 


If € 3n for real part, 0631 
No, evaluate f(z,,,) 
If é -n for imag part, 0725 
Build and 

store testword 
Reset A, pickup 

in O 
Load 


2141 
Load 


A 
Set exit in cpx. mpy. routine 
Transfer into cpx. mpy. 
Temp. store imag part 
Set exit in add routine 
Load A 
Transfer {nto add routine 
Load imag part for cpx. mpy. 
Reload 214] 

for next term 
Modify A, pickup 
If A, not taken, 

repeat back 
Testword storage 
Constant for 062 
Store f(24,1) 
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ADD 


0661 
0662 
0663 
0664 
0665 
0666 
0667 
0670 
0671 
0672 
0673 
0674 
0675 
0676 
0677 
0700 
0701 
0702 
0703 
0704 
0705 
0706 
0707 
0710 
0711 
0712 
0713 
0714 
0715 
0716 
0717 
0720 
0721 
0722 
0723 
0724 
0725 
0726 
0727 
0730 
0731 
0732 
0733 
0734 
0735 
0736 
0737 
0740 
0741 


INST 


M+ 


1307 
1327 
3000 
2000 
2002 
2100 
3000 
0676 
2000 
0677 
1330 
1332 
3000 
0000 
5200 
2010 
3000 
1310 
1351 
J 5g 
1255 
1501 
1521 
1502 
1522 
1510 
1512 
1550 
1332 
1306 
1326 
1307 
1327 
140 
1442 
3000 
1312 
0731 
0750 
0747 
0100 
2001 
2000 
2002 
0104 
3000 
2000 
2001 
0735 


M? 


1350 
1370 
2100 
2001 
2100 
2003 
2100 
2000 
0676 
2001 
1637 
1637 
2100 
0000 


3000 
2100 
2100 
2100 
2100 
2100 
1342 
1362 
1343 
1363 
1351 
15595 
1571 
1515 
1347 
1367 
1350 
1370 
1401 
1403 
2100 
1455 
1466 
1677 
1640 
0101 
1351 
1310 
2003 
0105 


2100 
2100 
1466 


M? 


2004 
2006 
1650 
2004 
2006 
2007 
1607 
0700 
0673 
0700 
1330 
1332 
0605 
0007 
0650 
0702 
0706 
0001 
0001 
0001 
0001 
1300 
1320 
1301 
1321 
1302 
1322 
1303 
1323 
1305 
1325 
1306 
1326 
1307 
1327 
0361 
0751 
0735 
0737 
0740 
2000 
2003 
2002 
2002 
2000 
1540 
0104 
0105 
0735 


IV-1-8 


Jana + 


REMARKS 


Load 
f(z, ) 

f (24,3 )/f (24) 

Load quotient 
for complex 
con Jugate 
multiplication 


Lf 
Ip (2444 )/f (2 |? 
is greater ar 
100; halve 

y and 
recompute 24, 

Floating point 
100 

If SEN 1 up, print iteration 

No, skip print 

Print 
214] 
real and 
imaginary 

Shift 21-1] 
to 24.0 

Shift 24 
to 21-1 

onift 2141 
to Zi 

Shift 44) 
to > 

Shift f(z ) 
to f(2z1. 

Shift f (24 
to 02 1 

Shift f(z1+1 ) 
to f(21) 

Reiterate 

If imag portion significant, 0751 

Initialize 0735 

Initialize 07357 

Initialize O740 

Load A, 


An 24 


An 21 
Put away 

terms of suppressed equation 
Modify to suppress 
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ADD 


07h2 
0743 
0744 
0745 
0746 
0747 
0750 
0751 
0752 
0753 
0754 
0755 
0756 
0757 
0760 
0761 
0762 
0763 
0764 
0765 
0766 
0767 
0770 
0771 
0772 
0773 
OTTh 
0775 
0776 
0777 
1000 
1001 
1002 
1003 
1004 
1005 
1006 
1007 
1010 
1011 
1012 
1013 
1014 
1015 
1016 
1017 
1020 
1021 


INST 


M1 


0737 
0740 
0656 
1372 
3000 
0000 
0000 
1310 
1351 
3000 
2002 
2000 
2100 
3000 
2000 
0766 
2007 
1457 
1002 
2007 
0100 
0102 
0114 
1300 
1341 
2000 
1300 
1341 
2000 
1301 
1302 
3000 
1300 
1341 
2000 
2002 
0770 
1002 
1003 
0656 
1372 
1351 
1310 
1551 
2100 
1351 
2004 
5000 


M? 


0264 
0264 
0735 
1331 
2100 
0000 
0000 
1551 
2100 
3000 
2005 
2100 
2001 
2100 
2001 
1473 
0264 
1677 
1331 
1953 
0101 
0103 
0115 
1303 
1344 
2001 
1304 
1345 
2001 
1342 
1343 
2100 
2100 
2100 
2001 
2003 
14.66 
0264 
0264 
0770 
0264 
1310 
0261 
2000 
2100 
1166 
2100 
2100 


M3 


OrsT 
0740 
0732 
1372 
1015 
0105 
0102 
1303 
1344 
1310 
2006 
2004 
2005 
1607 
1304 
2007 
0770 
1002 
2007 
1003 
1300 
1301 
1302 
2000 
2001 
2004 
2000 
2001 
2006 
2000 
2002 
1410 
0110 
0111 
1300 
1301 
0770 
1002 
100% 
0770 
15712 
1301 
2000 
2004 
1302 
1302 
1102 
1031 


REMARKS 


eqution to next 
lower degree 
If extraction not complete, 0752 
Reduce N' by 1 
Transfer to conversion routine 
Constant for 0730 
Constant for 0727 
Build 2a term 
in quadratic factor 
Build 
(a^ + 9%) 
term in 
quadratic 
factor 


Initialize 

0770 
Initialize 1002 
Initialize 

1003 
An to To 
Ans] to Fs b 


Ane? to po 


In x a to 2004, 2005 


A, x (a? + b?) to 2006, 2007 


Load 
An-] 
Add 
Putaway A, of 
suppressed equation 
Set up for 
next term 
Modify to pick up An.3 
Modify putaway 
Modify putaway 
Has extraction been completed? 
Yes, reduce N’ by 2 
Form fractional part (1) 
Build shift control word 
Form integer part 
Clear sign storage 
Store sign 
If integer > 0, 1102 
No, 1031 


(1) The remainder of Channel 10 and Channel 11 is a binary floating 
point to decimal conversion which would be part of the overall progran. 


1V = 1 - 9 








LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD 


1022 
1023 
1024 
1025 
1026 
1027 
1030 
1031 
1032 
1033 
1034 
1055 
1036 
1057 
1040 
1041 
1042 
1043 
1044 
1045 
1046 
1047 
1050 
1051 
1052 
1053 
1054 
1055 
1056 
1057 
1060 
1061 
1062 
1063 
1064 
1065 
1066 
1067 
1070 
1071 
1072 
1075 
1074 
1075 
1076 
1077 
1100 
1101 


INST 


12 
10 
00 
00 
10 
00 
30 
30 
30 
26 
27 
32 
37 
50 
32 
3l 
34 
34 
27 
32 
34 
21 
30 
30 
34 
32 
21 
21 
5h 
21 
32 
27 
34 
21 
34 
21 
21 
25 
21l 
59 
51 
3h 
21 
3l 
22 
55 
34 
01 


M+ 


0000 
0100 
2400 
1515 
0000 
3115 
2005 
1301 
1167 
2000 
2007 
2001 
2007 
2007 
1016 
1163 
3000 
1331 
1101 
1300 
2000 
1024 
1300 
2007 
3000 
1302 
1300 
1025 
3000 
1027 
1302 
2000 
2000 
1157 
3000 
1160 
1301 
1455 
1161 
1163 
1512 
3000 
1154 
1372 
0000 
1164 
3000 
0000 


Me 


0000 
0000 
0000 
2700 
0100 
1515 
2100 
2100 
2100 
0257 
0256 
0254 
2100 
2100 
1022 
1067 
2100 
2004 
2100 
0260 
2100 
1023 
0256 
0256 
2100 
1022 
2007 
1026 
2100 
1140 
1166 
0262 
2100 
1026 
2100 
1026 
1101 
1312 
1165 
2100 
T329 
2100 
1023 
2100 
0000 
2100 
2100 
0000 


IV 


M3 


0000 
0000 
0000 
0000 


1500 
1300 
2000 
2007 
2000 
2007 
2007 
1035 
1301 
1301 
1043 
1170 
1057 
2007 
2100 
1053 
0001 
1300 
2007 
1045 
1300 
0001 
0001 
1066 
0001 
2100 
2000 
1065 
0001 
1066 
0001 
0001 
1075 
0002 
1067 
1310 
1013 
0001 
0310 
0000 
1067 
1074 
0001 


- 1-10 


REMARKS 


Suppress sign bits 
Print control 
Alphabetic print 
Alphabetic print 
Print control 
Alphabetic print 
Reset command 
Load fraction 
Load testword 
Multiply for most sig. digit 
Shift testword 
Extract digit into TW storage 
Repeat loop if o'flo 
Putaway converted fraction 
Suppress sign 
Print editing 
commands 
Print editing 
commends 
Extract most sig. dec. digit 
If 2000 > 0, 1053 
If not print space 
Shift word hà left 
and change print control 
Return to test next digit 
Enter sign 
Print integer 
Print decimal point 
Skip to 1066 
Print editing cormand 
Build sign print 
when only fractional 
part involved 
Print +. 
Skip to 1066 
Print -. 
Print fractional part 
If imag part insig, 1074 
No, print + j 
Set skip command in 1067 
Set up convert 
imaginary part 
Print carriage return 
If N* 30, next root 
If N' = O, hait 
Restore 1067 
Return to "test for end" 
Constant for 104% 








LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD 


1102 
1103 
1104 
1105 
1106 
1107 
1110 
1111 
1112 
1113 
1114 
1115 
1116 
1117 
1120 
1121 
1122 
1123 
1124 
1125 
1126 
1127 
1130 
1131 
1132 
1133 
1134 
1135 
1136 
1137 
1140 
1141 
1142 
1143 
1144 
1145 
1146 
1147 
1150 
72351 
1152 
1153 
1154 
1155 
1156 


INST 


う 2 
31 
32 
30 
う 2 
35 
25 
26 
34 
02 
00 
00 
22 
25 
25 
23 
23 
23 
23 


M? 


2100 
1113 
2000 
2000 
1116 
1126 
2004 
2000 
3000 
0000 
0000 
0000 
1115 
2004 
2004 
2004 
2004 
2004 
2004 
2004 
2004 
0167 
0013 
0001 
0000 
0000 
0000 
3000 
0000 
0000 


2000 
2005 
1110 
2001 
1110 
3000 
0000 
0000 
2100 
3000 
2004 
5500 
2000 
0000 


M? 


2100 
2004 
111) 
1115 
2000 
2100 
1140 
1147 
2100 
0000 
0000 
0000 
2100 
1127 
1130 
1131 
1132 
1155 
1134 
1136 
1157 
1531 
2274 
1422 
0750 
0060 
0004 
2100 
0000 
0000 
0000 
1155 
1150 
1156 
2005 
11959 
2100 
0000 
0000 
2100 
2100 
1140 
0000 
0000 
0001 


M? 


2005 
2000 
2100 
2000 
1107 
1110 
2000 
2000 
1142 
0002 
0174 
0026 
1110 
2000 
2000 
2000 
2000 
2000 
2000 
2000 
2000 
LUCI 
0777 
2377 
e177 
6477 
1057 
1151 
0307 
0025 
0100 
2000 
2005 
1110 
2005 
1141 
1030 
0024 
000) 
2005 
1146 
2000 
0000 
0012 
0000 


IV =- 1 = 11 


High 


REMARKS 


integer 
conversion 
routine 








LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD INST mi M? M? REMARKS 

1157 00 2352 2764 6400 Print configuration for 1063 
1160 00 2252 2764 6400 Print configuration for 1065 
1161 00 2315 3621 3264 Print configuration for 1070 
1162 00 426% 6464 3232 Print configuration for 1070 
1163 34 3000 2100 1077 Dummy reset for 1071 

1164 33 1455 1312 1074 Durmy reset for 1074 

1165 10 0000 0000 0000 Print control - alphabetic 
1166 02 0000 0000 0000 Sign extractor 

1167 OO 0h21 0h21 0420 Testword for 1072 

1370 30 1016 1022 1302 Suppress integer sign when 
1171 324 2004 2100 1044 printing imaginary portion 
1172 21 1174 1175 0001 or print O. 

1175 5 3000 2100 1066 and return to 1066 

1174 00 5227 0000 0000 Print configuration for 1172 
1179 10 0001 0000 0000 Print control word 

1175 117, Not required 


Channel 12 is used for a memory check-sum routine in the root-solving 
routine and is available for other use when this program is incorporated 
in the overall synthesis program. 


Channel 13 is used for temporary storage and for a few constants which 
are listed below. NOTE: This program is coded for minimum access only, 
i.e., adjacent cells in the main memory are numbered 00, 41, 02, 43, 
O4, 45, etc. 


1551, 1557 - floating point "1", must be addressed individually 
1555, 1555 - floating point "O", must be addressed individually 


1572 - N' - the degree of the reduced equation 


1374 
1376 


00 


00 


3000 


0000 


3000 


0000 


0200 


0010 


Testword for 0305 


Increment for 0303, 0304 


Channels 14, 15, and 16 contain complex floating point arithmetic 
routines, 1.e., 


Complex Addition 1410 
Complex Subtraction 1416 
Simple Multiplication 1500 
Simple Division 1507 
Simple Subtraction 1516 
Simple Addition 1540 
Square Root 1551 
Complex Multiplication 1607 
Complex Division 1650 


IV =~ 1 = 12 








LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD INST 
1400 - 1507 
1410 30 
1111 32 
1412 Oh 
1413 32 
1414 32 
1415 34 
1416 30 
1417 32 
1420 0% 
1421 32 
1422 32 
1423 31 
1424 34 
1425 31 
1426 51 
1427 31 
1530 5h 
1431 31 
1432 31 
1433 34 
1434 = 1436 
1437 00 
1440 - 1453 
1454 02 
1455 - 1463 
1464 00 
1465 00 
1466 00 
1467 00 
1470 31 
1471 00 
1,72 00 
1h75 OO 
1474 00 
1475 - 1477 
1500 30 
1501 32 
1502 55 
1505 25 
1504 34 
1505 02 
1506 00 
1507 30 
1510 32 
1511 36 
1512 23 
15P5 37 
1514 34 
1515 00 


M! M? M? REMARKS 
Used for temporary storage of operands 
3000 1505 1400 Prepare 
1400 1640 1433 exit 
3000 3000 1400 Store operands 
1472 1535 1424 Set 1424 
1471 1677 1430 and 1430 for addition 
3000 2100 1423 Transfer to 1h25 
3000 1603 1536 Prepare 
1536 515 1433 exit 
3000 3000 1400 Store operands 
1465 1636 1424 Set 1424 
1467 1636 1430 and 1430 for subtraction 
2004 2005 2002 Reposition real operands 
3000 2100 1516 Transfer for + reals 
2000 2001 1400 Store sum of reals 
1402 1443 2000 Load imaginary 
1406 1447 2002 opefands 
3000 2100 1516 Transfer for + imaginaries 
2000 2001 2002 Reposition imaginary sum 
1400 1441 2000 Load real sum 
3000 2100 0624 Exit 
Temporary storage 
0000 0000 0100 Decrementer for 0306, 0307 
Temporary storage 
0000 0000 0044 
Temporary storage 
0000 7777 17777 M^, M? extractor for 0555 
0000 0000 1516 Constant for 1h21 
0002 0002 0000 Constant for 0323 
0000 0000 1516 Constant for 1422 
0102 0103 2000 Dummy reset for 0326 
0000 0000 1540 Constant for 1414 
0000 0000 1540 Constant for 1413 
0004 0004 0000 Constant for 0332, 0336 
0002 0000 0000 Incrementer (various) 
Temporary storage 
3000 1505 2006 Prepare 
2006 1677 1550 exit 
2000 2004 2000 Add exponents 
2001 2003 #2001 Multiply magnitudes 
3000 2100 1547 Transfer to exit 
0000 0000 0040 Shift control for 1500 
0000 0000 3777 M? extractor 
3000 1676 2006 Prepare 
2006 1677 1550 exit 
2000 2004 2000 Subtract exponents 
2001 2005 2001 Divide magnitudes 
2001 3000 1530 If overflow, 1530 
3000 2100 1547 Transfer to exot 
0000 0000 IE M extractor 


IV = 1 = 13 








LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 





ADD INST M! M? M? REMARKS 

1516 30 3000 1505 2006 Prepare 

1517 32 2006 1677 1550 exit 

1520 36 2100 2003 2003 Change sign 274 operand 

1521 34 3000 2100 1542 Transfer into add 

1522 36 2000 2002 2006 Form shift control 

1523 30 2001 2006 2007 Shift T operand 

1524 30 2002 2100 2000 Shift 2°° exponent to Ans. cell 

1525 55 2005 2007 2001 Add operands 

1526 $T 2001 3000 1530 If overflow, 1530 

1527 53h 3000 2100 1547 Transfer to exit 

1530 27 2001 1644 2001 Normal 

1531 30 2001 1644 2001 correct 

1532 27 2001 1646 2001 overflow 

1533 35 2000 1646 2000 procedure 

1534 3, 3000 2100 1547 Transfer to exit 

1535. 00 0000 0000 7777 M? extractor 

1536 00 0000 0000 06214 Prepare exit storage 

1537 02 0000 0000 0001 Shift control for 1555 

1540 30 3000 1505 2006 Prepare 

1541 32 2006 1515 1550 exit 

1542 33 2002 2000 1522 If 2002 > 2000, 1522 

1543 36 2002 2000 2006 No, form, shift control 

1544 30 2003 2006 2007 Shift 224 operand 

1545 35 2001 2007 2001 Add operands 

1546 37 2001 3000 1530 If overflow, 1530 

1547 31l 2000 2001 2000 Normalize 

1550 34 3000 2100 E Exit 

1551 50 3000 1676 2006 Prepare exit 

1552 . $2 2006 1955 1550 from square root routine 

1553 35 200 0261 2000 Multiply by 2** 

1554 32 2000 7551 2104 Extract last bit of exp into 200k 

1555 30 2000 1644 2006 Take square root of exp 

1556 34 2004 2100 1561 If exp odd, magnitude has been 
divided by 2, by 1555 

1557 30 2001 1644 2007 If exp even, divide magnitude 
by 2 

1560 3 3000 2100 1562 Skip to iteration 

1561 40 2001 2100 2007 load magnitude for iteration 

1562 40 1575 2100 2002 Load ay 

1563 23 2007 2002 2001 (x/2)/a, 

1564 36 2001 2002 2001 (x/2a y - 8, 

1565 30 2001 164 2001 1/2(x/2a,) = - 84 = 2441 - B4 

1566 34 2001 2100 1570 Has process converged? 

1567 34 3000 2100 1572 Yes, 1572; no, 1570 

1570 35 2002 2001 2002 aj] replaces aj 

1571 5h 3000 2100 1565 and reiterate 

1572 25 2002 1577 2001 Multilpy by 42/2 

1573 30 2006 1576 2000 Divide by 2?! 

1574 34 3000 2100 1547 Transfer to exit 


IV -1-1h 








LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 


ADD 


1379 
1576 
D 
1600 
1603 
1604 
1607 
1610 
1611 
1612 
1613 
1614 
1615 
1616 
1617 
1620 
1621 
1622 
1623 
1624 
1625 
1626 
1627 
1630 
1631 
1632 
1633 
1634 
1635 
1636 
1637 
1640 
1641 
1642 
1643 
1644 
1645 
1646 
1647 
1650 
1651 
1652 
1653 
1654 
1655 


Ine x E 


INST Mi M? M? 
00 Te TT T 
00 0000 0000 0021 
00 5520 2363 1500 
1602 Used 
02 0000 0000 0030 
1606 Used 
30 3000 1603 1601 
32 1601 1636 1635 
ok 3000 3000 1600 
34 3000 2100 1500 
30 2000 2100 1601 
30 2001 2100 1605 
35 1602 1606 2000 
25 1643 1647 2001 
31 2000 2001 2002 
21 1601 1605 2000 
34 3000 2100 1516 
30 2000 2100 1601 
30 2001 2100 1605 
35 1602 1604 2000 
25 1643 1645 2001 
31 2000 2001 2002 
う 5 1600 1606 2000 
25 1641 1647 2001 
3] 2000 2001 2000 
34 3000 2100 1540 
3] 2000 2001 2002 
21 1601 1605 2000 
3l 3000 2100 0627 
00 0000 0000 NT 
OO 0000 0000 0001 
00 0000 0000 TEE 
Used for temporary storage 
00 0000 0000 0001 
Used for temporary storage 
02 0000 0000 0001 
Used for temporary storage 
00 0000 0000 0001 
Used for temporary storage 
30 3000 1676 1601 
32 1601 1677 1635 
0l 3000 3000 1600 
30 2004 1637 2000 
25 2005 1645 2001 
31 2000 2001 2000 


REMARKS 


Constant for 1562 
Constant for 1573 
Constant for 1572 


for temporary storage of operands 


Shift control for 1607 


for temporary storage of operands 


Prepare 
exit 
Store operands 
Multiply reals 
and store results 
in 1601, 1605 
Multiply imeginaries 
to obtain negative real 
and shift for subtraction 
from previous real product 
Transfer to subtract 
Putaway real product 
in 1601, 1605 
Multiply for 
first cross-product term 
and hold in 2002, 2003 
Multiply second 
cross-product term 
and prepare for addition 
Add for imaginary cross-products 
Put results in2002, 2003 
Insert real product 
Exit 
M^ extractor for 1610 
Shifter for 1656 
extractor 


Shift control word 

Shift control word 

Shift Control word 

Prepare (A +4B) 
exit (C +3D) 

Store operands 

"orm 


c? 
end normalire 
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ADD 


1656 
1657 
1660 
1661 
1662 
1663 
1664 
1665 
1666 
1667 
1670 
1671 
1572 
1673 
157% 
1675 
1676 
1677 


INST 


20 
25 
3] 
34 
31 
3] 
34 
30 
30 
31 
34 
30 
36 
3] 


mM? 


2006 
2007 
2002 
3000 
2000 
1604 
3000 
2000 
2001 
1606 
3000 
2000 
2100 
1604 
1600 
3000 
0000 
0000 


M? 


1637 
1647 
2003 
2100 
2001 
1645 
2100 
2100 
2100 
1647 
2100 
2100 
2001 
1645 
1641 
2100 
0000 
0000 


FINIS 


IV - 1-16 


M? 


2002 
2003 
2002 
1540 
2004 
2000 
1507 
1604 
1645 
2000 
1507 
1606 
1647 
2004 
2000 
1612 
0030 
DIT 


REMARKS 
Form 


and normalize 
(Ca DF) 
Shift to 2004, 2005 
Load C 
CINC ) 
Store in 
1604, 1645 
Load D 
pus. pe) 
Store in 1606, 
1647 with negative sign 
Load conjugate 
quotient and 
transfer to multiply 
Shift control word 
extractor 
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APPENDIX V 


USE OF THE EXISTING PROGRAMS ON THE CRC-102A 


All existing programs have been coded for minimum-access word channel 
and will not work in the sequential channels. 
l. Dirichlet Power Series Expansion 

Since the synthetic division and conversion program has not been 
written as yet, if £ (t) and f (t) are given, Int A must be computed by 
hand. If h(t) is given, it must be converted to In a by multiplying by 
At. As the program is now located, storage is available for 100 (octal) 
ordinates in floating point form if the decimal to floating point con- 
version in Channel 2 is to be used. If the conversion is accomplished 


prior to entry of data, 150 (octal) ordinates may be used. 


S = 2 — me 2NT where S = total cells needed 


T = number of time units problem is 
scaled over (note command 0423) 


At = sampling interval 
N = samples per time unit 

The program may be relocated to accommodate more sampled ordinates. 

The binary floating point numbers representing the ordinates of h(t) 
are entered in cells 0000 - O11/7, exponents in even cells, magnitudes in 
odd; this allows 40 ordinates to be used. 

The test switches determine the amount of interpolation performed. 
The output will be H*(s) in decimal form for as many terms requested. 
The number of terms is set by entering this datum as an octal number in 
J, cell 0540. 
2. Matrix Inversion 

The matrix inversion program takes the ascending power H*(s) terms in 


Vol 





packed floating point form (a conversion routine is available) in the word 
structure shown on page III-5 and computes the lossless system function, 
Z'(p). The terms are entered in Channel 0,0000, 0001, 0002, etc. Cells 
0076 and 0077 are reserved for N, the order of the matrix which equals n, 
the number of poles desired or the degree of the denominator and n + r, 
the degree of the numerator respectively. For the example on page 21, 
Nene4 (0076) 
n*rzsn-22 (0077) 
The commands 0300 through 0320 contain instructions which link the program 
to a packed binary floating point to decimal conversion routine located in 
1571 - 1677. 

If these commands are included, the program will yield the decimal 
coefficients of Z'p in the following form: 

Po + PS + pos” 

do + 915 + q,5^ - ds + EN 
do will be equal to 1, the results must be normalized to the following 
form for entry in the root-solving routine. 

(M) a. + pis * p. 
l 0 
E + TS + a + q!s +q? 
3 2 1 0 

3. Lagrangian Interpolation Method 

The procedure for the use of this program is clearly described in 
the description of sub-routines for the Computation Center (Routine 2.1). 

It is recommended that the Channel 2 conversion routine of this pro- 
gram be replaced to provide the normalizing procedure required by 2. above. 
Channels 11 and 12 are available for the storage of the factored forms of 


the numerator and denominator. 
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